diff --git a/Framework/PythonInterface/plugins/algorithms/CalculateEfficiencyCorrection.py b/Framework/PythonInterface/plugins/algorithms/CalculateEfficiencyCorrection.py index ba81c727c9d0c2ee1535b1bd19c3a75d6b818b0e..5020f5004a8f5c40728fa09ebc0b77876473c076 100644 --- a/Framework/PythonInterface/plugins/algorithms/CalculateEfficiencyCorrection.py +++ b/Framework/PythonInterface/plugins/algorithms/CalculateEfficiencyCorrection.py @@ -8,17 +8,19 @@ from __future__ import absolute_import, division, print_function import numpy as np from mantid.simpleapi import \ - CloneWorkspace, ConvertFromDistribution, ConvertToPointData, ConvertUnits, SetSampleMaterial + CloneWorkspace, ConvertFromDistribution, ConvertToPointData, \ + ConvertUnits, CreateWorkspace, Rebin, SetSampleMaterial from mantid.api import \ - mtd, AlgorithmFactory, DataProcessorAlgorithm, MatrixWorkspaceProperty + AlgorithmFactory, CommonBinsValidator, PropertyMode, PythonAlgorithm, \ + MatrixWorkspaceProperty from mantid.kernel import \ - Direction, FloatBoundedValidator, MaterialBuilder, StringMandatoryValidator + Direction, FloatBoundedValidator, MaterialBuilder, \ + StringListValidator, StringMandatoryValidator -PROPS_FOR_TRANS = ['DensityType'] TABULATED_WAVELENGTH = 1.7982 -class CalculateEfficiencyCorrection(DataProcessorAlgorithm): +class CalculateEfficiencyCorrection(PythonAlgorithm): _input_ws = None _output_ws = None @@ -34,30 +36,48 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): or vanadium measurements.' def seeAlso(self): - return [ "He3TubeEfficiency", "DetectorEfficiencyCor", "DetectorEfficiencyCorUser", + return [ "He3TubeEfficiency", "CalculateSampleTransmission", + "DetectorEfficiencyCor", "DetectorEfficiencyCorUser", "CalculateEfficiency", "ComputeCalibrationCoefVan" ] def PyInit(self): self.declareProperty( MatrixWorkspaceProperty('InputWorkspace', '', - direction=Direction.Input), + direction=Direction.Input, + optional=PropertyMode.Optional, + validator=CommonBinsValidator()), doc='Input workspace with wavelength range to calculate for the correction.') + self.declareProperty(name='WavelengthRange', defaultValue='', + doc='Wavelength range to calculate efficiency for.') + self.declareProperty( MatrixWorkspaceProperty('OutputWorkspace', '', direction=Direction.Output), doc="Outputs the efficiency correction which can be applied by multiplying \ the OutputWorkspace to the workspace that requires the correction.") + self.declareProperty( + name='ChemicalFormula', defaultValue='None', + validator=StringMandatoryValidator(), + doc='Sample chemical formula used to determine :math:`\\sigma (\\lambda_{ref})` term.') + + self.declareProperty( + name='DensityType', defaultValue = 'Mass Density', + validator=StringListValidator(['Mass Density', 'Number Density']), + doc = 'Use of Mass density or Number density') + self.declareProperty( name='Density', defaultValue=0.0, + validator=FloatBoundedValidator(0.0), doc='Mass density (g/cm^3) or Number density (atoms/Angstrom^3), :math:`\\rho`. \ Default=0.0') self.declareProperty( name='Thickness', defaultValue=1.0, + validator=FloatBoundedValidator(0.0), doc='Sample thickness (cm), :math:`T`. Default value=1.0') self.declareProperty( @@ -65,9 +85,13 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): defaultValue=0.0, validator=FloatBoundedValidator(0.0, 1.0), doc="Directly input the efficiency, :math:`\\epsilon`, measured at \ - MeasuredEfficiencyWavelength, :math:`\\lambda_{\\epsilon}`, to determine \ - :math:`\\rho * T = - ln(1-\\epsilon) \\frac{\\lambda_{ref}}{\\lambda_{\\epsilon} \\sigma_{a}(\\lambda_{ref})}` term, \ - where :math:`\\lambda_{ref} =` 1.7982. Default=0.0") + MeasuredEfficiencyWavelength, :math:`\\lambda_{\\epsilon}`, to determine :math:`\\rho * T`, where \ + :math:`\\rho * T = - ln(1-\\epsilon) \\frac{1}{ \\frac{\\lambda_{\\epsilon} \\sigma (\\lambda_{ref})}{\\lambda_{ref}}}` \ + term for XSectionType == AttenuationXSection and \ + :math:`\\rho * T = - ln(1-\\epsilon) \\frac{1} \ + { \\sigma_s + \\frac{\\lambda_{\\epsilon} \\sigma (\\lambda_{ref})}{\\lambda_{ref}}}` \ + for XSectionType == TotalXSection \ + with :math:`\\lambda_{ref} =` 1.7982. Default=0.0") self.declareProperty( name='MeasuredEfficiencyWavelength', @@ -80,25 +104,27 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): name='Alpha', defaultValue=0.0, doc="Directly input :math:`\\alpha`, where \ - :math:`\\alpha = \\rho * T * \\sigma_{a}(\\lambda_{ref}) / \\lambda_{ref}`, \ - where :math:`\\lambda_{ref} =` 1.7982. Default=0.0") + :math:`\\alpha = \\rho * T * \\sigma (\\lambda_{ref}) / \\lambda_{ref}`, \ + where :math:`\\lambda_{ref} =` 1.7982. XSectionType has no effect. Default=0.0") self.declareProperty( - name='ChemicalFormula', defaultValue='None', - validator=StringMandatoryValidator(), - doc='Sample chemical formula used to determine :math:`\\sigma_{a}(\\lambda_{ref}) term.`') - - self.copyProperties('CalculateSampleTransmission', PROPS_FOR_TRANS) + name='XSectionType', + defaultValue="AttenuationXSection", + validator=StringListValidator(['AttenuationXSection', 'TotalXSection']), + doc = 'Use either the absorption cross section (for monitor-type correction) or, \ + the total cross section (for transmission-type correction) in Alpha term. \ + Default=AttenuationXSection') def validateInputs(self): issues = dict() - density = self.getProperty('Density').value - if density < 0.0: - issues['Density'] = 'Density must be positive' - thickness = self.getProperty('Thickness').value - if thickness < 0.0: - issues['Thickness'] = 'Thickness must be positive' + if self.getProperty('InputWorkspace').isDefault and self.getProperty('WavelengthRange').isDefault: + issues['InputWorkspace'] = "Must select either InputWorkspace and WavelengthRange as input" + issues['WavelengthRange'] = "Must select either InputWorkspace and WavelengthRange as input" + + if not self.getProperty('InputWorkspace').isDefault and not self.getProperty('WavelengthRange').isDefault: + issues['InputWorkspace'] = "Cannot select both InputWorkspace and WavelengthRange as input" + issues['WavelengthRange'] = "Cannot select both InputWorkspace and WavelengthRange as input" if not self.getProperty('Density').isDefault: if self.getProperty('ChemicalFormula').isDefault: @@ -112,15 +138,15 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): if self.getProperty('ChemicalFormula').isDefault: issues['ChemicalFormula'] = "Must specify the ChemicalFormula with MeasuredEfficiency" - if not self.getProperty('MeasuredEfficiency').isDefault and not self.getProperty('Density'): + if not self.getProperty('MeasuredEfficiency').isDefault and not self.getProperty('Density').isDefault: issues['MeasuredEfficiency'] = "Cannot select both MeasuredEfficiency and Density as input" issues['Density'] = "Cannot select both MeasuredEfficiency and Density as input" - if not self.getProperty('Alpha').isDefault and not self.getProperty('Density'): + if not self.getProperty('Alpha').isDefault and not self.getProperty('Density').isDefault: issues['Alpha'] = "Cannot select both Alpha and Density as input" issues['Density'] = "Cannot select both Alpha and Density as input" - if not self.getProperty('MeasuredEfficiency').isDefault and not self.getProperty('Alpha'): + if not self.getProperty('MeasuredEfficiency').isDefault and not self.getProperty('Alpha').isDefault: issues['MeasuredEfficiency'] = "Cannot select both MeasuredEfficiency and Alpha as input" issues['Alpha'] = "Cannot select both MeasuredEfficiency and Alpha as input" @@ -128,90 +154,102 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): def _setup(self): self._input_ws = self.getPropertyValue('InputWorkspace') - self._efficiency = self.getProperty('MeasuredEfficiency').value - self._efficiency_wavelength = self.getProperty('MeasuredEfficiencyWavelength').value + self._bin_params = self.getPropertyValue('WavelengthRange') + self._output_ws = self.getProperty('OutputWorkspace').valueAsStr self._chemical_formula = self.getPropertyValue('ChemicalFormula') self._density_type = self.getPropertyValue('DensityType') self._density = self.getProperty('Density').value self._thickness = self.getProperty('Thickness').value - self._alpha = self.getProperty('Alpha').value - self._output_ws = self.getProperty('OutputWorkspace').valueAsStr + self._efficiency = self.getProperty('MeasuredEfficiency').value + self._efficiency_wavelength = self.getProperty('MeasuredEfficiencyWavelength').value + self._alpha_absXS = self.getProperty('Alpha').value + self._alpha_scatXS = 0.0 + self._xsection_type = self.getProperty('XSectionType').value def PyExec(self): self._setup() - CloneWorkspace( - Inputworkspace=self._input_ws, - OutputWorkspace=self._output_ws) - - if mtd[self._output_ws].isDistribution(): - ConvertFromDistribution(Workspace=self._output_ws) - - ConvertToPointData( - InputWorkspace=self._output_ws, - OutputWorkspace=self._output_ws) - ConvertUnits( - InputWorkspace=self._output_ws, - OutputWorkspace=self._output_ws, - Target='Wavelength', - EMode='Elastic') + if not self.getProperty("InputWorkspace").isDefault: + self._output_ws = CloneWorkspace(Inputworkspace=self._input_ws, + OutputWorkspace=self._output_ws, + StoreInADS=False) + else: + self._output_ws = CreateWorkspace(NSpec=1, DataX=[0], DataY=[0], + UnitX='Wavelength', Distribution=False, + StoreInADS=False) + self._output_ws = Rebin(InputWorkspace=self._output_ws, + Params=self._bin_params, + StoreInADS=False) + + if self._output_ws.isDistribution(): + ConvertFromDistribution(Workspace=self._output_ws, + StoreInADS=False) + + self._output_ws = ConvertToPointData(InputWorkspace=self._output_ws, + StoreInADS=False) + self._output_ws = ConvertUnits(InputWorkspace=self._output_ws, + Target='Wavelength', + EMode='Elastic', + StoreInADS=False) if self.getProperty('Alpha').isDefault: SetSampleMaterial( InputWorkspace=self._output_ws, ChemicalFormula=self._chemical_formula, - SampleNumberDensity=self._density) + SampleNumberDensity=self._density, + StoreInADS=False) if self.getProperty('MeasuredEfficiency').isDefault: self._calculate_area_density_from_density() else: self._calculate_area_density_from_efficiency() - self._calculate_alpha() + self._calculate_alpha_absXS_term() + if self._xsection_type == "TotalXSection": + self._calculate_alpha_scatXS_term() - ws = mtd[self._output_ws] - wavelengths = ws.readX(0) + wavelengths = self._output_ws.readX(0) efficiency = self._calculate_efficiency(wavelengths) - for histo in range(ws.getNumberHistograms()): - mtd[self._output_ws].setY(histo, efficiency) + for histo in range(self._output_ws.getNumberHistograms()): + self._output_ws.setY(histo, efficiency) - self.setProperty('OutputWorkspace', mtd[self._output_ws]) + self.setProperty('OutputWorkspace', self._output_ws) def _calculate_area_density_from_efficiency(self): - ''' - Calculates area density (atom/cm^2) using efficiency - - ''' - material = mtd[self._output_ws].sample().getMaterial() + """Calculates area density (atom/cm^2) using efficiency""" + material = self._output_ws.sample().getMaterial() ref_absXS = material.absorbXSection() - self._area_density = - np.log( 1.0 - self._efficiency) / ref_absXS - self._area_density *= TABULATED_WAVELENGTH / self._efficiency_wavelength + xs_term = ref_absXS * self._efficiency_wavelength / TABULATED_WAVELENGTH + if self._xsection_type == "TotalXSection": + xs_term += material.totalScatterXSection() + self._area_density = - np.log( 1.0 - self._efficiency) / xs_term def _calculate_area_density_from_density(self): - ''' - Calculates area density (atom/cm^2) using number density and thickness. - - ''' + """Calculates area density (atom/cm^2) using number density and thickness.""" if self._density_type == 'Mass Density': builder = MaterialBuilder().setFormula(self._chemical_formula) mat = builder.setMassDensity(self._density).build() self._density = mat.numberDensity self._area_density = self._density * self._thickness - def _calculate_alpha(self): - ''' - Calculates alpha = -area_density / absXS / 1.7982 - ''' + def _calculate_alpha_absXS_term(self): + """Calculates absorption XS alpha term = area_density * absXS / 1.7982""" + material = self._output_ws.sample().getMaterial() + absXS = material.absorbXSection() / TABULATED_WAVELENGTH + self._alpha_absXS = self._area_density * absXS - material = mtd[self._output_ws].sample().getMaterial() - absorption_x_section = material.absorbXSection() / TABULATED_WAVELENGTH - self._alpha = self._area_density * absorption_x_section + def _calculate_alpha_scatXS_term(self): + """Calculates scattering XS alpha term = area_density * scatXS""" + material = self._output_ws.sample().getMaterial() + scatXS = material.totalScatterXSection() + self._alpha_scatXS = self._area_density * scatXS def _calculate_efficiency(self, wavelength): - ''' + """ Calculates efficiency of a detector / monitor at a given wavelength. + If using just the absorption cross section, _alpha_scatXS is == 0. @param wavelength Wavelength at which to calculate (in Angstroms) @return efficiency - ''' - efficiency = 1. / (1.0 - np.exp(-self._alpha * wavelength)) + """ + efficiency = 1. / (1.0 - np.exp(-(self._alpha_scatXS + self._alpha_absXS * wavelength))) return efficiency diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/CalculateEfficiencyCorrectionTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/CalculateEfficiencyCorrectionTest.py index 0a838b16bfd5278e4a4a05fe7a985f8b6e657499..503bcdd479746587f82d5fe423b00e3fdf6f905b 100644 --- a/Framework/PythonInterface/test/python/plugins/algorithms/CalculateEfficiencyCorrectionTest.py +++ b/Framework/PythonInterface/test/python/plugins/algorithms/CalculateEfficiencyCorrectionTest.py @@ -8,7 +8,8 @@ from __future__ import (absolute_import, division, print_function) import unittest from mantid.simpleapi import \ - LoadAscii, DeleteWorkspace, Multiply, ConvertToPointData, CalculateEfficiencyCorrection + ConvertToPointData, CalculateEfficiencyCorrection, CreateSampleWorkspace, \ + DeleteWorkspace, LoadAscii, Multiply from testhelpers import run_algorithm from mantid.api import AnalysisDataService @@ -18,12 +19,15 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): _input_wksp = "input_wksp" _correction_wksp = "correction_wksp" _output_wksp = "output_wksp" + _wavelengths="0.2,0.01,4.0" _alpha = 0.693 _chemical_formula = "(He3)" _number_density = 0.0002336682167635477 _mass_density = 0.0011702649036052439 - _efficiency1 = 0.712390781371 - _efficiency2 = { "Efficiency": 0.74110947758, "Wavelength": 1.95} + _efficiency1_forAbsXS = 0.712390781371 + _efficiency2_forAbsXS = { "Efficiency": 0.74110947758, "Wavelength": 1.95} + _efficiency1_forTotalXS = 0.712793729636 + _efficiency2_forTotalXS = { "Efficiency": 0.741472190178, "Wavelength": 1.95} _thickness = 1.0 def setUp(self): @@ -39,13 +43,12 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): S. Howells, "On the choice of moderator for a liquids diffractometer on a pulsed neutron source", Nucl. Instr. Meth. Phys. Res. 223, 1984, doi:10.1016/0167-5087(84)90256-4 ''' - input_wksp = LoadAscii( - Filename="CalculateEfficiencyCorrection_milder_moderator_polyethlyene_300K.txt", - Unit="Wavelength") - ConvertToPointData(InputWorkspace=input_wksp, OutputWorkspace=input_wksp) - self._input_wksp = input_wksp + LoadAscii(OutputWorkspace=self._input_wksp, + Filename="CalculateEfficiencyCorrection_milder_moderator_polyethlyene_300K.txt", + Unit="Wavelength") + ConvertToPointData(InputWorkspace=self._input_wksp, OutputWorkspace=self._input_wksp) - def checkResults(self): + def checkResults(self, xsection="AttenuationXSection"): Multiply(LHSWorkspace=self._input_wksp, RHSWorkspace=self._correction_wksp, OutputWorkspace=self._output_wksp) @@ -53,9 +56,43 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): self.assertEqual(output_wksp.getAxis(0).getUnit().unitID(), 'Wavelength') self.assertAlmostEqual(output_wksp.readX(0)[79], 0.995) - self.assertAlmostEqual(output_wksp.readY(0)[79], 3250.28183501) + if xsection == "AttenuationXSection": + self.assertAlmostEqual(output_wksp.readY(0)[79], 3250.28183501) + if xsection == "TotalXSection": + self.assertAlmostEqual(output_wksp.readY(0)[79], 3245.70148939) # Result tests + def testCalculateEfficiencyCorrectionAlphaForEventWksp(self): + self._input_wksp = "input_wksp" + self._correction_wksp = "correction_wksp" + self._output_wksp = "output_wksp" + + # Create an exponentially decaying function in wavelength to simulate + # measured sample + CreateSampleWorkspace(WorkspaceType="Event", Function="User Defined", + UserDefinedFunction="name=ExpDecay, Height=100, Lifetime=4", + Xmin=0.2, Xmax=4.0, BinWidth=0.01, XUnit="Wavelength", + NumEvents=10000, NumBanks=1, BankPixelWidth=1, + OutputWorkspace=self._input_wksp) + + # Calculate the efficiency correction + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + Alpha=self._alpha, + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + # Apply correction + ConvertToPointData(InputWorkspace=self._input_wksp, OutputWorkspace=self._input_wksp) + Multiply(LHSWorkspace=self._input_wksp, + RHSWorkspace=self._correction_wksp, + OutputWorkspace=self._output_wksp) + + # Check results + output_wksp = AnalysisDataService.retrieve(self._output_wksp) + self.assertEqual(output_wksp.getAxis(0).getUnit().unitID(), 'Wavelength') + self.assertAlmostEqual(output_wksp.readX(0)[79], 0.995) + self.assertAlmostEqual(output_wksp.readY(0)[79], 62.22517501) def testCalculateEfficiencyCorrectionAlpha(self): alg_test = run_algorithm("CalculateEfficiencyCorrection", @@ -66,63 +103,259 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): self.checkResults() + def testCalculateEfficiencyCorrectionAlphaWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + Alpha=self._alpha, + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + def testCalculateEfficiencyCorrectionMassDensityAndThickness(self): - correction_wksp = "correction_ws" alg_test = run_algorithm("CalculateEfficiencyCorrection", InputWorkspace=self._input_wksp, ChemicalFormula=self._chemical_formula, Density=self._mass_density, Thickness=self._thickness, - OutputWorkspace=correction_wksp) + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionMassDensityAndThicknessWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + Density=self._mass_density, + Thickness=self._thickness, + OutputWorkspace=self._correction_wksp) self.assertTrue(alg_test.isExecuted()) - self._correction_wksp = AnalysisDataService.retrieve(correction_wksp) self.checkResults() def testCalculateEfficiencyCorrectionNumberDensityAndThickness(self): - correction_wksp = "correction_ws" alg_test = run_algorithm("CalculateEfficiencyCorrection", InputWorkspace=self._input_wksp, ChemicalFormula=self._chemical_formula, DensityType="Number Density", Density=self._number_density, Thickness=self._thickness, - OutputWorkspace=correction_wksp) + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionNumberDensityAndThicknessWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + DensityType="Number Density", + Density=self._number_density, + Thickness=self._thickness, + OutputWorkspace=self._correction_wksp) self.assertTrue(alg_test.isExecuted()) - self._correction_wksp = AnalysisDataService.retrieve(correction_wksp) self.checkResults() def testCalculateEfficiencyCorrectionEfficiency(self): - correction_wksp = "correction_ws" alg_test = run_algorithm("CalculateEfficiencyCorrection", InputWorkspace=self._input_wksp, ChemicalFormula=self._chemical_formula, - MeasuredEfficiency=self._efficiency1, - OutputWorkspace=correction_wksp) + MeasuredEfficiency=self._efficiency1_forAbsXS, + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionEfficiencyWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=self._efficiency1_forAbsXS, + OutputWorkspace=self._correction_wksp) self.assertTrue(alg_test.isExecuted()) - self._correction_wksp = AnalysisDataService.retrieve(correction_wksp) self.checkResults() def testCalculateEfficiencyCorrectionEfficiencyWithWavelength(self): - correction_wksp = "correction_ws" - efficiency=self._efficiency2['Efficiency'] - wavelength=self._efficiency2['Wavelength'] + efficiency=self._efficiency2_forAbsXS['Efficiency'] + wavelength=self._efficiency2_forAbsXS['Wavelength'] alg_test = run_algorithm("CalculateEfficiencyCorrection", InputWorkspace=self._input_wksp, ChemicalFormula=self._chemical_formula, MeasuredEfficiency=efficiency, MeasuredEfficiencyWavelength=wavelength, - OutputWorkspace=correction_wksp) + OutputWorkspace=self._correction_wksp) self.assertTrue(alg_test.isExecuted()) - self._correction_wksp = AnalysisDataService.retrieve(correction_wksp) self.checkResults() + def testCalculateEfficiencyCorrectionEfficiencyWithWavelengthWithWaveRange(self): + efficiency=self._efficiency2_forAbsXS['Efficiency'] + wavelength=self._efficiency2_forAbsXS['Wavelength'] + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=efficiency, + MeasuredEfficiencyWavelength=wavelength, + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionAlphaWithTotalXS(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + Alpha=self._alpha, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionAlphaWithTotalXSWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + Alpha=self._alpha, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults() + + def testCalculateEfficiencyCorrectionMassDensityAndThicknessWithTotalXS(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + ChemicalFormula=self._chemical_formula, + Density=self._mass_density, + Thickness=self._thickness, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionMassDensityAndThicknessWithTotalXSWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + Density=self._mass_density, + Thickness=self._thickness, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionNumberDensityAndThicknessWithTotalXS(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + ChemicalFormula=self._chemical_formula, + DensityType="Number Density", + Density=self._number_density, + Thickness=self._thickness, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionNumberDensityAndThicknessWithTotalXSWithWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + DensityType="Number Density", + Density=self._number_density, + Thickness=self._thickness, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionEfficiencyWithTotalXS(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=self._efficiency1_forTotalXS, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionEfficiencyWithTotalXSWitWaveRange(self): + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=self._efficiency1_forTotalXS, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionEfficiencyWithWavelengthWithTotalXS(self): + efficiency=self._efficiency2_forTotalXS['Efficiency'] + wavelength=self._efficiency2_forTotalXS['Wavelength'] + alg_test = run_algorithm("CalculateEfficiencyCorrection", + InputWorkspace=self._input_wksp, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=efficiency, + MeasuredEfficiencyWavelength=wavelength, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorrectionEfficiencyWithWavelengthWithTotalXSWithWaveRange(self): + efficiency=self._efficiency2_forTotalXS['Efficiency'] + wavelength=self._efficiency2_forTotalXS['Wavelength'] + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + ChemicalFormula=self._chemical_formula, + MeasuredEfficiency=efficiency, + MeasuredEfficiencyWavelength=wavelength, + XSectionType="TotalXSection", + OutputWorkspace=self._correction_wksp) + self.assertTrue(alg_test.isExecuted()) + + self.checkResults(xsection="TotalXSection") + + def testCalculateEfficiencyCorretionStoreADSCheck(self): + self.cleanup() + alg_test = run_algorithm("CalculateEfficiencyCorrection", + WavelengthRange=self._wavelengths, + Alpha=self._alpha, + OutputWorkspace=self._output_wksp) + self.assertTrue(alg_test.isExecuted()) + self.assertTrue(AnalysisDataService.doesExist(self._output_wksp)) + # Invalid checks - def testCalculateEfficiencyCorretionInvalidDensity(self): + def testCalculateEfficiencyCorretionInvalidStoreADSCheck(self): + self.cleanup() + corr_wksp = CalculateEfficiencyCorrection( + WavelengthRange=self._wavelengths, + Alpha=self._alpha, + StoreInADS=False) + self.assertFalse(AnalysisDataService.doesExist(corr_wksp.name())) + + def testCalculateEfficiencyCorretionInvalidInput(self): self.assertRaises(RuntimeError, + CalculateEfficiencyCorrection, + OutputWorkspace=self._output_wksp) + + def testCalculateEfficiencyCorretionInvalidTooManyInputs(self): + self.assertRaises(RuntimeError, + CalculateEfficiencyCorrection, + InputWorkspace=self._input_wksp, + WavelengthRange=self._wavelengths, + OutputWorkspace=self._output_wksp) + + def testCalculateEfficiencyCorretionInvalidDensity(self): + self.assertRaises(ValueError, CalculateEfficiencyCorrection, InputWorkspace=self._input_wksp, ChemicalFormula=self._chemical_formula, @@ -175,11 +408,13 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): MeasuredEfficiency=1.0, OutputWorkspace=self._output_wksp) - def cleanUp(self): + def cleanup(self): if AnalysisDataService.doesExist(self._input_wksp): DeleteWorkspace(self._input_wksp) if AnalysisDataService.doesExist(self._output_wksp): DeleteWorkspace(self._output_wksp) + if AnalysisDataService.doesExist(self._correction_wksp): + DeleteWorkspace(self._correction_wksp) if __name__ == "__main__": diff --git a/docs/source/algorithms/CalculateEfficiencyCorrection-v1.rst b/docs/source/algorithms/CalculateEfficiencyCorrection-v1.rst index f3e934c197a47ae4e3bd7ca475026a4800848a08..6db84f168b380bdbec3240358a726bb5e89ac449 100644 --- a/docs/source/algorithms/CalculateEfficiencyCorrection-v1.rst +++ b/docs/source/algorithms/CalculateEfficiencyCorrection-v1.rst @@ -14,8 +14,8 @@ This algorithm calculates the detection efficiency correction, :math:`\epsilon`, .. math:: :label: efficiency - \epsilon &= 1-e^{-\rho T \sigma_a (\lambda)} \\ - &= 1-e^{-\rho_{A} \sigma_a (\lambda)} \\ + \epsilon &= 1-e^{-\rho T \sigma (\lambda)} \\ + &= 1-e^{-\rho_{A} \sigma (\lambda)} \\ &= 1-e^{-\alpha \lambda} and output correction is given as: @@ -25,25 +25,32 @@ and output correction is given as: \frac{1}{\epsilon} = \frac{1}{1-e^{-\alpha \lambda}} -where :math:`\rho` is the number density in :math:`atoms/\AA^3`, -:math:`T` is a sample thickness given in :math:`cm`, -:math:`\sigma_a(\lambda)` is the wavelength-dependent absorption cross-section given as -:math:`\sigma_a (\lambda) = \sigma_a (\lambda_{ref}) \left( \frac{\lambda}{\lambda_{ref}} \right)`, -in units of :math:`barns` and with :math:`\lambda_{ref}` = 1.7982 :math:`\AA`, -:math:`\rho_{A}` is the area density (:math:`\rho_{A}=\rho * T`) in units of :math:`atoms*cm/\AA^3`, -:math:`\alpha = \rho_{A} * \frac{\sigma_a (\lambda_{ref})}{\lambda_{ref}} = \rho * T * \frac{\sigma_a (\lambda_{ref})}{\lambda_{ref}}` in units of :math:`1/\AA`. -and :math:`\lambda` is in units of :math:`\AA`. +where -NOTE: :math:`1 \AA^2 = 10^{8} barns` and :math:`1 \AA = 10^{-8} cm`. +- :math:`\rho` is the number density in atoms/:math:`\AA^3` +- :math:`T` is a sample thickness given in cm +- :math:`\lambda_{ref}` = 1.7982 :math:`\AA`, +- :math:`\sigma (\lambda)` is the wavelength-dependent cross-section which is either: + + - :math:`\sigma (\lambda) = \sigma_a (\lambda_{ref}) \left( \frac{\lambda}{\lambda_{ref}} \right)` for ``XSectionType`` == ``AttenuationXSection`` where :math:`\sigma_a` is the absorption cross-section in units of barns + - or :math:`\sigma (\lambda) = \sigma_s + \sigma_a (\lambda_{ref}) \left( \frac{\lambda}{\lambda_{ref}} \right)` for ``XSectionType`` == ``TotalXSection`` where :math:`\sigma_s` is the total scattering cross-section in units of barns + +- :math:`\rho_{A}` is the area density (:math:`\rho_{A}=\rho * T`) in units of atoms*cm/:math:`\AA^3`, +- :math:`\alpha = \rho_{A} * \frac{\sigma (\lambda_{ref})}{\lambda_{ref}} = \rho * T * \frac{\sigma (\lambda_{ref})}{\lambda_{ref}}` in units of 1/:math:`\AA`. +- :math:`\lambda` is in units of :math:`\AA`. + +NOTE: :math:`1 \AA^2 = 10^{8}` barns and :math:`1 \AA = 10^{-8}` cm. The optional inputs into the algorithm are to input either: 1. The ``Alpha`` parameter - 2. The ``Density`` and ``ChemicalFormula`` to calculate the :math:`\sigma_a(\lambda_{ref})` term. - 3. The ``MeasuredEfficiency``, an optional ``MeasuredEfficiencyWavelength``, and ``ChemicalFormula`` to calculate the :math:`\sigma_a(\lambda_{ref})` term. + 2. The ``Density`` and ``ChemicalFormula`` to calculate the :math:`\sigma(\lambda_{ref})` term. + 3. The ``MeasuredEfficiency``, an optional ``MeasuredEfficiencyWavelength``, and ``ChemicalFormula`` to calculate the :math:`\sigma(\lambda_{ref})` term. The ``MeasuredEfficiency`` is the :math:`\epsilon` term measured at a specific wavelength, ``MeasuredEfficiencyWavelength``. This helps if the efficiency has been directly measured experimentally at a given wavelength. +For the ``XSectionType``, if the efficiency correction is applied to a beam monitor to determine the incident spectrum, then the ``AttenuationXSection`` option should be used. This is due to the fact that scatter events do not lead to neutrons that will be in the incident beam. If the efficiency correction is to be used similar to a transmission measurement for an actual sample measurement, such as in :ref:`algm-CalculateSampleTransmission-v1`, then the ``TotalXSection`` should be used to include both types of events. + Usage ----- **Example - Basics of running CalculateEfficiencyCorrection with Alpha.** @@ -57,15 +64,18 @@ Usage NumEvents=10000, NumBanks=1, BankPixelWidth=1) # Calculate the efficiency correction - correction_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, Alpha=0.5) + corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, Alpha=0.5) + corr_wksp_with_wave_range = CalculateEfficiencyCorrection(WavelengthRange="0.2,0.01,4.0", Alpha=0.5) # Apply the efficiency correction to the measured spectrum input_wksp = ConvertToPointData(InputWorkspace=input_wksp) - output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=correction_wksp) + output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp) + output_wksp_with_wave_range = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp_with_wave_range) print('Input workspace: {}'.format(input_wksp.readY(0)[:5])) - print('Correction workspace: {}'.format(correction_wksp.readY(0)[:5])) + print('Correction workspace: {}'.format(corr_wksp.readY(0)[:5])) print('Output workspace: {}'.format(output_wksp.readY(0)[:5])) + print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5])) Ouptut: @@ -74,6 +84,7 @@ Ouptut: Input workspace: [ 38. 38. 38. 38. 38.] Correction workspace: [ 10.26463773 9.81128219 9.39826191 9.02042771 8.67347109] Output workspace: [ 390.05623383 372.82872321 357.13395265 342.77625306 329.59190131] + Output workspace using WavelengthRange: [ 390.05623383 372.82872321 357.13395265 342.77625306 329.59190131] **Example - Basics of running CalculateEfficiencyCorrection with Density and ChemicalFormula.** @@ -86,17 +97,22 @@ Ouptut: NumEvents=10000, NumBanks=1, BankPixelWidth=1) # Calculate the efficiency correction - correction_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, - Density=6.11, - ChemicalFormula="V") + corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, + Density=6.11, + ChemicalFormula="V") + corr_wksp_with_wave_range = CalculateEfficiencyCorrection(WavelengthRange="0.2,0.01,4.0", + Density=6.11, + ChemicalFormula="V") # Apply the efficiency correction to the measured spectrum input_wksp = ConvertToPointData(InputWorkspace=input_wksp) - output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=correction_wksp) + output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp) + output_wksp_with_wave_range = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp_with_wave_range) print('Input workspace: {}'.format(input_wksp.readY(0)[:5])) - print('Correction workspace: {}'.format(correction_wksp.readY(0)[:5])) + print('Correction workspace: {}'.format(corr_wksp.readY(0)[:5])) print('Output workspace: {}'.format(output_wksp.readY(0)[:5])) + print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5])) Ouptut: @@ -105,6 +121,7 @@ Ouptut: Input workspace: [ 38. 38. 38. 38. 38.] Correction workspace: [ 24.40910309 23.29738394 22.28449939 21.35783225 20.50682528] Output workspace: [ 927.54591732 885.30058981 846.81097679 811.59762534 779.25936055] + Output workspace using WavelengthRange: [ 927.54591732 885.30058981 846.81097679 811.59762534 779.25936055] **Example - Basics of running CalculateEfficiencyCorrection with MeasuredEfficiency and ChemicalFormula.** @@ -117,17 +134,24 @@ Ouptut: NumEvents=10000, NumBanks=1, BankPixelWidth=1) # Calculate the efficiency correction - correction_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, - MeasuredEfficiency=1e-2, - ChemicalFormula="(He3)") + corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, + MeasuredEfficiency=1e-2, + ChemicalFormula="(He3)") + + corr_wksp_with_wave_range = CalculateEfficiencyCorrection(WavelengthRange="0.2,0.01,4.0", + MeasuredEfficiency=1e-2, + ChemicalFormula="(He3)") + # Apply the efficiency correction to the measured spectrum input_wksp = ConvertToPointData(InputWorkspace=input_wksp) - output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=correction_wksp) + output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp) + output_wksp_with_wave_range = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp_with_wave_range) print('Input workspace: {}'.format(input_wksp.readY(0)[:5])) - print('Correction workspace: {}'.format(correction_wksp.readY(0)[:5])) + print('Correction workspace: {}'.format(corr_wksp.readY(0)[:5])) print('Output workspace: {}'.format(output_wksp.readY(0)[:5])) + print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5])) Ouptut: @@ -137,6 +161,83 @@ Ouptut: Correction workspace: [ 873.27762699 832.68332786 795.69741128 761.85923269 730.78335476] Output workspace: [ 33184.54982567 31641.9664586 30236.50162877 28950.65084207 27769.76748099] + Output workspace using WavelengthRange: [ 33184.54982567 31641.9664586 30236.50162877 28950.65084207 + 27769.76748099] + +**Example - Basics of running CalculateEfficiencyCorrection with MeasuredEfficiency and ChemicalFormula using the total cross section.** + +.. testcode:: ExBasicCalcualteEfficiencyCorrectionWithEfficiency + + # Create an exponentially decaying function in wavelength to simulate measured sample + input_wksp = CreateSampleWorkspace(WorkspaceType="Event", Function="User Defined", + UserDefinedFunction="name=ExpDecay, Height=100, Lifetime=4", + Xmin=0.2, Xmax=4.0, BinWidth=0.01, XUnit="Wavelength", + NumEvents=10000, NumBanks=1, BankPixelWidth=1) + + # Calculate the efficiency correction + corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, + MeasuredEfficiency=1e-2, + ChemicalFormula="(He3)", + XSectionType="TotalXSection") + + corr_wksp_with_wave_range = CalculateEfficiencyCorrection(WavelengthRange="0.2,0.01,4.0", + MeasuredEfficiency=1e-2, + ChemicalFormula="(He3)", + XSectionType="TotalXSection") + + + # Apply the efficiency correction to the measured spectrum + input_wksp = ConvertToPointData(InputWorkspace=input_wksp) + output_wksp = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp) + output_wksp_with_wave_range = Multiply(LHSWorkspace=input_wksp, RHSWorkspace=corr_wksp_with_wave_range) + + print('Input workspace: {}'.format(input_wksp.readY(0)[:5])) + print('Correction workspace: {}'.format(corr_wksp.readY(0)[:5])) + print('Output workspace: {}'.format(output_wksp.readY(0)[:5])) + print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5])) + +Ouptut: + +.. testoutput:: ExBasicCalcualteEfficiencyCorrectionWithEfficiency + + Input workspace: [ 38. 38. 38. 38. 38.] + Correction workspace: [ 865.7208838 825.85320701 789.49774383 756.20995361 725.61727932] + Output workspace: [ 32897.39358441 31382.42186624 30000.91426562 28735.97823706 + 27573.45661411] + Output workspace using WavelengthRange: [ 32897.39358441 31382.42186624 30000.91426562 28735.97823706 + 27573.45661411] + +The transmission of a sample can be measured as :math:`e^{-\rho T \sigma_t (\lambda)}` where :math:`\sigma_t (\lambda) = \sigma_s + \sigma_a (\lambda)` is the total cross-section. This can be calculatd directly by the :ref:`algm-CalculateSampleTransmission-v1` algorithm. Yet, we can also back out the transmission with the ``CalculateEfficiencyCorrection`` algorithm. The example below shows how: + +**Example - Transmission using the CalculateEfficiencyCorrection and CalculateSampleTransmission comparison.** + +.. testcode:: ExTransmissionCalcualteEfficiencyCorrection + + ws = CalculateSampleTransmission(WavelengthRange='2.0, 0.1, 10.0', + ChemicalFormula='H2-O') + print('Transmission: {} ...'.format(ws.readY(0)[:3])) + + corr_wksp = CalculateEfficiencyCorrection(WavelengthRange="2.0, 0.1, 10.0", + Density=0.1, + Thickness=0.1, + ChemicalFormula="H2-O", + XSectionType="TotalXSection") + dataX = corr_wksp.readX(0) + dataY = np.ones(len(corr_wksp.readX(0))) + ones = CreateWorkspace(dataX, dataY, UnitX="Wavelength") + efficiency = Divide(LHSWorkspace=ones, RHSWorkspace=corr_wksp) # 1 + -1 * transmission + negative_trans = Minus(LHSWorkspace=efficiency, RHSWorkspace=ones) # -1 * transmission + transmission = Multiply(LHSWOrkspace=negative_trans, RHSWorkspace=-1.*ones) + print('Transmission using efficiency correction: {} ...'.format(transmission.readY(0)[:3])) + +Output: + +.. testoutput:: ExTransmissionCalcualteEfficiencyCorrection + + Transmission: [ 0.94506317 0.94505148 0.94503979] ... + Transmission using efficiency correction: [ 0.9450632 0.94505151 0.94503982] ... + +The discrepancies are due to the differenc in :math:`\lambda_{ref}` = 1.7982 :math:`\AA` in ``CalculateEfficiencyCorrection``, consistent with `ReferenceLambda <https://github.com/mantidproject/mantid/blob/32ed0b2cbbe4fbfb230570d5a53032f6101743de/Framework/Kernel/src/NeutronAtom.cpp#L23>`_ where :ref:`algm-CalculateSampleTransmission-v1` uses :math:`\lambda_{ref}` = 1.798 :math:`\AA`. **Example - Running CalculateEfficiencyCorrection for incident spectrum.**