Skip to content
Snippets Groups Projects
Commit 3caec64b authored by McDonnell, Marshall's avatar McDonnell, Marshall
Browse files

Refactor for CalculateEfficiencyCorrection algorithm Refs #24432

parent b09e102f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/PythonInterface/plugins/algorithms/CalculateEfficiencyCorrection.py
+ 111
− 73
View file @ 3caec64b
...@@ -8,17 +8,19 @@ from __future__ import absolute_import, division, print_function ...@@ -8,17 +8,19 @@ from __future__ import absolute_import, division, print_function
import numpy as np import numpy as np
from mantid.simpleapi import \ from mantid.simpleapi import \
CloneWorkspace, ConvertFromDistribution, ConvertToPointData, ConvertUnits, SetSampleMaterial CloneWorkspace, ConvertFromDistribution, ConvertToPointData, \
ConvertUnits, CreateWorkspace, Rebin, SetSampleMaterial
from mantid.api import \ from mantid.api import \
mtd, AlgorithmFactory, DataProcessorAlgorithm, MatrixWorkspaceProperty AlgorithmFactory, CommonBinsValidator, PropertyMode, PythonAlgorithm, \
MatrixWorkspaceProperty
from mantid.kernel import \ from mantid.kernel import \
Direction, FloatBoundedValidator, MaterialBuilder, StringMandatoryValidator Direction, FloatBoundedValidator, MaterialBuilder, \
StringListValidator, StringMandatoryValidator
PROPS_FOR_TRANS = ['DensityType']
TABULATED_WAVELENGTH = 1.7982 TABULATED_WAVELENGTH = 1.7982
class CalculateEfficiencyCorrection(DataProcessorAlgorithm): class CalculateEfficiencyCorrection(PythonAlgorithm):
_input_ws = None _input_ws = None
_output_ws = None _output_ws = None
...@@ -34,30 +36,48 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): ...@@ -34,30 +36,48 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm):
or vanadium measurements.' or vanadium measurements.'
def seeAlso(self): def seeAlso(self):
return [ "He3TubeEfficiency", "DetectorEfficiencyCor", "DetectorEfficiencyCorUser", return [ "He3TubeEfficiency", "CalculateSampleTransmission",
"DetectorEfficiencyCor", "DetectorEfficiencyCorUser",
"CalculateEfficiency", "ComputeCalibrationCoefVan" ] "CalculateEfficiency", "ComputeCalibrationCoefVan" ]
def PyInit(self): def PyInit(self):
self.declareProperty( self.declareProperty(
MatrixWorkspaceProperty('InputWorkspace', '', MatrixWorkspaceProperty('InputWorkspace', '',
direction=Direction.Input), direction=Direction.Input,
optional=PropertyMode.Optional,
validator=CommonBinsValidator()),
doc='Input workspace with wavelength range to calculate for the correction.') 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( self.declareProperty(
MatrixWorkspaceProperty('OutputWorkspace', '', MatrixWorkspaceProperty('OutputWorkspace', '',
direction=Direction.Output), direction=Direction.Output),
doc="Outputs the efficiency correction which can be applied by multiplying \ doc="Outputs the efficiency correction which can be applied by multiplying \
the OutputWorkspace to the workspace that requires the correction.") 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( self.declareProperty(
name='Density', name='Density',
defaultValue=0.0, defaultValue=0.0,
validator=FloatBoundedValidator(0.0),
doc='Mass density (g/cm^3) or Number density (atoms/Angstrom^3), :math:`\\rho`. \ doc='Mass density (g/cm^3) or Number density (atoms/Angstrom^3), :math:`\\rho`. \
Default=0.0') Default=0.0')
self.declareProperty( self.declareProperty(
name='Thickness', name='Thickness',
defaultValue=1.0, defaultValue=1.0,
validator=FloatBoundedValidator(0.0),
doc='Sample thickness (cm), :math:`T`. Default value=1.0') doc='Sample thickness (cm), :math:`T`. Default value=1.0')
self.declareProperty( self.declareProperty(
...@@ -65,9 +85,13 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): ...@@ -65,9 +85,13 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm):
defaultValue=0.0, defaultValue=0.0,
validator=FloatBoundedValidator(0.0, 1.0), validator=FloatBoundedValidator(0.0, 1.0),
doc="Directly input the efficiency, :math:`\\epsilon`, measured at \ doc="Directly input the efficiency, :math:`\\epsilon`, measured at \
MeasuredEfficiencyWavelength, :math:`\\lambda_{\\epsilon}`, to determine \ MeasuredEfficiencyWavelength, :math:`\\lambda_{\\epsilon}`, to determine :math:`\\rho * T`, where \
:math:`\\rho * T = - ln(1-\\epsilon) \\frac{\\lambda_{ref}}{\\lambda_{\\epsilon} \\sigma_{a}(\\lambda_{ref})}` term, \ :math:`\\rho * T = - ln(1-\\epsilon) \\frac{1}{ \\frac{\\lambda_{\\epsilon} \\sigma (\\lambda_{ref})}{\\lambda_{ref}}}` \
where :math:`\\lambda_{ref} =` 1.7982. Default=0.0") 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( self.declareProperty(
name='MeasuredEfficiencyWavelength', name='MeasuredEfficiencyWavelength',
...@@ -80,25 +104,27 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): ...@@ -80,25 +104,27 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm):
name='Alpha', name='Alpha',
defaultValue=0.0, defaultValue=0.0,
doc="Directly input :math:`\\alpha`, where \ doc="Directly input :math:`\\alpha`, where \
:math:`\\alpha = \\rho * T * \\sigma_{a}(\\lambda_{ref}) / \\lambda_{ref}`, \ :math:`\\alpha = \\rho * T * \\sigma (\\lambda_{ref}) / \\lambda_{ref}`, \
where :math:`\\lambda_{ref} =` 1.7982. Default=0.0") where :math:`\\lambda_{ref} =` 1.7982. XSectionType has no effect. Default=0.0")
self.declareProperty( self.declareProperty(
name='ChemicalFormula', defaultValue='None', name='XSectionType',
validator=StringMandatoryValidator(), defaultValue="AttenuationXSection",
doc='Sample chemical formula used to determine :math:`\\sigma_{a}(\\lambda_{ref}) term.`') validator=StringListValidator(['AttenuationXSection', 'TotalXSection']),
doc = 'Use either the absorption cross section (for monitor-type correction) or, \
self.copyProperties('CalculateSampleTransmission', PROPS_FOR_TRANS) the total cross section (for transmission-type correction) in Alpha term. \
Default=AttenuationXSection')
def validateInputs(self): def validateInputs(self):
issues = dict() issues = dict()
density = self.getProperty('Density').value
if density < 0.0:
issues['Density'] = 'Density must be positive'
thickness = self.getProperty('Thickness').value if self.getProperty('InputWorkspace').isDefault and self.getProperty('WavelengthRange').isDefault:
if thickness < 0.0: issues['InputWorkspace'] = "Must select either InputWorkspace and WavelengthRange as input"
issues['Thickness'] = 'Thickness must be positive' 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 not self.getProperty('Density').isDefault:
if self.getProperty('ChemicalFormula').isDefault: if self.getProperty('ChemicalFormula').isDefault:
...@@ -112,15 +138,15 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm): ...@@ -112,15 +138,15 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm):
if self.getProperty('ChemicalFormula').isDefault: if self.getProperty('ChemicalFormula').isDefault:
issues['ChemicalFormula'] = "Must specify the ChemicalFormula with MeasuredEfficiency" 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['MeasuredEfficiency'] = "Cannot select both MeasuredEfficiency and Density as input"
issues['Density'] = "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['Alpha'] = "Cannot select both Alpha and Density as input"
issues['Density'] = "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['MeasuredEfficiency'] = "Cannot select both MeasuredEfficiency and Alpha as input"
issues['Alpha'] = "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): ...@@ -128,90 +154,102 @@ class CalculateEfficiencyCorrection(DataProcessorAlgorithm):
def _setup(self): def _setup(self):
self._input_ws = self.getPropertyValue('InputWorkspace') self._input_ws = self.getPropertyValue('InputWorkspace')
self._efficiency = self.getProperty('MeasuredEfficiency').value self._bin_params = self.getPropertyValue('WavelengthRange')
self._efficiency_wavelength = self.getProperty('MeasuredEfficiencyWavelength').value self._output_ws = self.getProperty('OutputWorkspace').valueAsStr
self._chemical_formula = self.getPropertyValue('ChemicalFormula') self._chemical_formula = self.getPropertyValue('ChemicalFormula')
self._density_type = self.getPropertyValue('DensityType') self._density_type = self.getPropertyValue('DensityType')
self._density = self.getProperty('Density').value self._density = self.getProperty('Density').value
self._thickness = self.getProperty('Thickness').value self._thickness = self.getProperty('Thickness').value
self._alpha = self.getProperty('Alpha').value self._efficiency = self.getProperty('MeasuredEfficiency').value
self._output_ws = self.getProperty('OutputWorkspace').valueAsStr 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): def PyExec(self):
self._setup() self._setup()
CloneWorkspace( if not self.getProperty("InputWorkspace").isDefault:
Inputworkspace=self._input_ws, self._output_ws = CloneWorkspace(Inputworkspace=self._input_ws,
OutputWorkspace=self._output_ws) OutputWorkspace=self._output_ws,
StoreInADS=False)
if mtd[self._output_ws].isDistribution(): else:
ConvertFromDistribution(Workspace=self._output_ws) self._output_ws = CreateWorkspace(NSpec=1, DataX=[0], DataY=[0],
UnitX='Wavelength', Distribution=False,
ConvertToPointData( StoreInADS=False)
InputWorkspace=self._output_ws, self._output_ws = Rebin(InputWorkspace=self._output_ws,
OutputWorkspace=self._output_ws) Params=self._bin_params,
ConvertUnits( StoreInADS=False)
InputWorkspace=self._output_ws,
OutputWorkspace=self._output_ws, if self._output_ws.isDistribution():
Target='Wavelength', ConvertFromDistribution(Workspace=self._output_ws,
EMode='Elastic') 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: if self.getProperty('Alpha').isDefault:
SetSampleMaterial( SetSampleMaterial(
InputWorkspace=self._output_ws, InputWorkspace=self._output_ws,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
SampleNumberDensity=self._density) SampleNumberDensity=self._density,
StoreInADS=False)
if self.getProperty('MeasuredEfficiency').isDefault: if self.getProperty('MeasuredEfficiency').isDefault:
self._calculate_area_density_from_density() self._calculate_area_density_from_density()
else: else:
self._calculate_area_density_from_efficiency() 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 = self._output_ws.readX(0)
wavelengths = ws.readX(0)
efficiency = self._calculate_efficiency(wavelengths) efficiency = self._calculate_efficiency(wavelengths)
for histo in range(ws.getNumberHistograms()): for histo in range(self._output_ws.getNumberHistograms()):
mtd[self._output_ws].setY(histo, efficiency) 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): def _calculate_area_density_from_efficiency(self):
''' """Calculates area density (atom/cm^2) using efficiency"""
Calculates area density (atom/cm^2) using efficiency material = self._output_ws.sample().getMaterial()
'''
material = mtd[self._output_ws].sample().getMaterial()
ref_absXS = material.absorbXSection() ref_absXS = material.absorbXSection()
self._area_density = - np.log( 1.0 - self._efficiency) / ref_absXS xs_term = ref_absXS * self._efficiency_wavelength / TABULATED_WAVELENGTH
self._area_density *= TABULATED_WAVELENGTH / self._efficiency_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): 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': if self._density_type == 'Mass Density':
builder = MaterialBuilder().setFormula(self._chemical_formula) builder = MaterialBuilder().setFormula(self._chemical_formula)
mat = builder.setMassDensity(self._density).build() mat = builder.setMassDensity(self._density).build()
self._density = mat.numberDensity self._density = mat.numberDensity
self._area_density = self._density * self._thickness self._area_density = self._density * self._thickness
def _calculate_alpha(self): def _calculate_alpha_absXS_term(self):
''' """Calculates absorption XS alpha term = area_density * absXS / 1.7982"""
Calculates alpha = -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() def _calculate_alpha_scatXS_term(self):
absorption_x_section = material.absorbXSection() / TABULATED_WAVELENGTH """Calculates scattering XS alpha term = area_density * scatXS"""
self._alpha = self._area_density * absorption_x_section material = self._output_ws.sample().getMaterial()
scatXS = material.totalScatterXSection()
self._alpha_scatXS = self._area_density * scatXS
def _calculate_efficiency(self, wavelength): def _calculate_efficiency(self, wavelength):
''' """
Calculates efficiency of a detector / monitor at a given 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) @param wavelength Wavelength at which to calculate (in Angstroms)
@return efficiency @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 return efficiency
......
Framework/PythonInterface/test/python/plugins/algorithms/CalculateEfficiencyCorrectionTest.py
+ 262
− 27
View file @ 3caec64b
...@@ -8,7 +8,8 @@ from __future__ import (absolute_import, division, print_function) ...@@ -8,7 +8,8 @@ from __future__ import (absolute_import, division, print_function)
import unittest import unittest
from mantid.simpleapi import \ from mantid.simpleapi import \
LoadAscii, DeleteWorkspace, Multiply, ConvertToPointData, CalculateEfficiencyCorrection ConvertToPointData, CalculateEfficiencyCorrection, CreateSampleWorkspace, \
DeleteWorkspace, LoadAscii, Multiply
from testhelpers import run_algorithm from testhelpers import run_algorithm
from mantid.api import AnalysisDataService from mantid.api import AnalysisDataService
...@@ -18,12 +19,15 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): ...@@ -18,12 +19,15 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase):
_input_wksp = "input_wksp" _input_wksp = "input_wksp"
_correction_wksp = "correction_wksp" _correction_wksp = "correction_wksp"
_output_wksp = "output_wksp" _output_wksp = "output_wksp"
_wavelengths="0.2,0.01,4.0"
_alpha = 0.693 _alpha = 0.693
_chemical_formula = "(He3)" _chemical_formula = "(He3)"
_number_density = 0.0002336682167635477 _number_density = 0.0002336682167635477
_mass_density = 0.0011702649036052439 _mass_density = 0.0011702649036052439
_efficiency1 = 0.712390781371 _efficiency1_forAbsXS = 0.712390781371
_efficiency2 = { "Efficiency": 0.74110947758, "Wavelength": 1.95} _efficiency2_forAbsXS = { "Efficiency": 0.74110947758, "Wavelength": 1.95}
_efficiency1_forTotalXS = 0.712793729636
_efficiency2_forTotalXS = { "Efficiency": 0.741472190178, "Wavelength": 1.95}
_thickness = 1.0 _thickness = 1.0
def setUp(self): def setUp(self):
...@@ -39,13 +43,12 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): ...@@ -39,13 +43,12 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase):
S. Howells, "On the choice of moderator for a liquids diffractometer on a pulsed neutron source", 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 Nucl. Instr. Meth. Phys. Res. 223, 1984, doi:10.1016/0167-5087(84)90256-4
''' '''
input_wksp = LoadAscii( LoadAscii(OutputWorkspace=self._input_wksp,
Filename="CalculateEfficiencyCorrection_milder_moderator_polyethlyene_300K.txt", Filename="CalculateEfficiencyCorrection_milder_moderator_polyethlyene_300K.txt",
Unit="Wavelength") Unit="Wavelength")
ConvertToPointData(InputWorkspace=input_wksp, OutputWorkspace=input_wksp) ConvertToPointData(InputWorkspace=self._input_wksp, OutputWorkspace=self._input_wksp)
self._input_wksp = input_wksp
def checkResults(self): def checkResults(self, xsection="AttenuationXSection"):
Multiply(LHSWorkspace=self._input_wksp, Multiply(LHSWorkspace=self._input_wksp,
RHSWorkspace=self._correction_wksp, RHSWorkspace=self._correction_wksp,
OutputWorkspace=self._output_wksp) OutputWorkspace=self._output_wksp)
...@@ -53,9 +56,43 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): ...@@ -53,9 +56,43 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase):
self.assertEqual(output_wksp.getAxis(0).getUnit().unitID(), 'Wavelength') self.assertEqual(output_wksp.getAxis(0).getUnit().unitID(), 'Wavelength')
self.assertAlmostEqual(output_wksp.readX(0)[79], 0.995) 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 # 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): def testCalculateEfficiencyCorrectionAlpha(self):
alg_test = run_algorithm("CalculateEfficiencyCorrection", alg_test = run_algorithm("CalculateEfficiencyCorrection",
...@@ -66,63 +103,259 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): ...@@ -66,63 +103,259 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase):
self.checkResults() 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): def testCalculateEfficiencyCorrectionMassDensityAndThickness(self):
correction_wksp = "correction_ws"
alg_test = run_algorithm("CalculateEfficiencyCorrection", alg_test = run_algorithm("CalculateEfficiencyCorrection",
InputWorkspace=self._input_wksp, InputWorkspace=self._input_wksp,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
Density=self._mass_density, Density=self._mass_density,
Thickness=self._thickness, 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.assertTrue(alg_test.isExecuted())
self._correction_wksp = AnalysisDataService.retrieve(correction_wksp)
self.checkResults() self.checkResults()
def testCalculateEfficiencyCorrectionNumberDensityAndThickness(self): def testCalculateEfficiencyCorrectionNumberDensityAndThickness(self):
correction_wksp = "correction_ws"
alg_test = run_algorithm("CalculateEfficiencyCorrection", alg_test = run_algorithm("CalculateEfficiencyCorrection",
InputWorkspace=self._input_wksp, InputWorkspace=self._input_wksp,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
DensityType="Number Density", DensityType="Number Density",
Density=self._number_density, Density=self._number_density,
Thickness=self._thickness, 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.assertTrue(alg_test.isExecuted())
self._correction_wksp = AnalysisDataService.retrieve(correction_wksp)
self.checkResults() self.checkResults()
def testCalculateEfficiencyCorrectionEfficiency(self): def testCalculateEfficiencyCorrectionEfficiency(self):
correction_wksp = "correction_ws"
alg_test = run_algorithm("CalculateEfficiencyCorrection", alg_test = run_algorithm("CalculateEfficiencyCorrection",
InputWorkspace=self._input_wksp, InputWorkspace=self._input_wksp,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
MeasuredEfficiency=self._efficiency1, MeasuredEfficiency=self._efficiency1_forAbsXS,
OutputWorkspace=correction_wksp) 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.assertTrue(alg_test.isExecuted())
self._correction_wksp = AnalysisDataService.retrieve(correction_wksp)
self.checkResults() self.checkResults()
def testCalculateEfficiencyCorrectionEfficiencyWithWavelength(self): def testCalculateEfficiencyCorrectionEfficiencyWithWavelength(self):
correction_wksp = "correction_ws" efficiency=self._efficiency2_forAbsXS['Efficiency']
efficiency=self._efficiency2['Efficiency'] wavelength=self._efficiency2_forAbsXS['Wavelength']
wavelength=self._efficiency2['Wavelength']
alg_test = run_algorithm("CalculateEfficiencyCorrection", alg_test = run_algorithm("CalculateEfficiencyCorrection",
InputWorkspace=self._input_wksp, InputWorkspace=self._input_wksp,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
MeasuredEfficiency=efficiency, MeasuredEfficiency=efficiency,
MeasuredEfficiencyWavelength=wavelength, MeasuredEfficiencyWavelength=wavelength,
OutputWorkspace=correction_wksp) OutputWorkspace=self._correction_wksp)
self.assertTrue(alg_test.isExecuted()) self.assertTrue(alg_test.isExecuted())
self._correction_wksp = AnalysisDataService.retrieve(correction_wksp)
self.checkResults() 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 # 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, 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, CalculateEfficiencyCorrection,
InputWorkspace=self._input_wksp, InputWorkspace=self._input_wksp,
ChemicalFormula=self._chemical_formula, ChemicalFormula=self._chemical_formula,
...@@ -175,11 +408,13 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase): ...@@ -175,11 +408,13 @@ class CalculateEfficiencyCorrectionTest(unittest.TestCase):
MeasuredEfficiency=1.0, MeasuredEfficiency=1.0,
OutputWorkspace=self._output_wksp) OutputWorkspace=self._output_wksp)
def cleanUp(self): def cleanup(self):
if AnalysisDataService.doesExist(self._input_wksp): if AnalysisDataService.doesExist(self._input_wksp):
DeleteWorkspace(self._input_wksp) DeleteWorkspace(self._input_wksp)
if AnalysisDataService.doesExist(self._output_wksp): if AnalysisDataService.doesExist(self._output_wksp):
DeleteWorkspace(self._output_wksp) DeleteWorkspace(self._output_wksp)
if AnalysisDataService.doesExist(self._correction_wksp):
DeleteWorkspace(self._correction_wksp)
if __name__ == "__main__": if __name__ == "__main__":
......
docs/source/algorithms/CalculateEfficiencyCorrection-v1.rst
+ 127
− 26
View file @ 3caec64b
...@@ -14,8 +14,8 @@ This algorithm calculates the detection efficiency correction, :math:`\epsilon`, ...@@ -14,8 +14,8 @@ This algorithm calculates the detection efficiency correction, :math:`\epsilon`,
.. math:: .. math::
:label: efficiency :label: efficiency
\epsilon &= 1-e^{-\rho T \sigma_a (\lambda)} \\ \epsilon &= 1-e^{-\rho T \sigma (\lambda)} \\
&= 1-e^{-\rho_{A} \sigma_a (\lambda)} \\ &= 1-e^{-\rho_{A} \sigma (\lambda)} \\
&= 1-e^{-\alpha \lambda} &= 1-e^{-\alpha \lambda}
and output correction is given as: and output correction is given as:
...@@ -25,25 +25,32 @@ 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}} \frac{1}{\epsilon} = \frac{1}{1-e^{-\alpha \lambda}}
where :math:`\rho` is the number density in :math:`atoms/\AA^3`, where
: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`.
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: The optional inputs into the algorithm are to input either:
1. The ``Alpha`` parameter 1. The ``Alpha`` parameter
2. The ``Density`` 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_a(\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 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. 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 Usage
----- -----
**Example - Basics of running CalculateEfficiencyCorrection with Alpha.** **Example - Basics of running CalculateEfficiencyCorrection with Alpha.**
...@@ -57,15 +64,18 @@ Usage ...@@ -57,15 +64,18 @@ Usage
NumEvents=10000, NumBanks=1, BankPixelWidth=1) NumEvents=10000, NumBanks=1, BankPixelWidth=1)
# Calculate the efficiency correction # 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 # Apply the efficiency correction to the measured spectrum
input_wksp = ConvertToPointData(InputWorkspace=input_wksp) 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('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: {}'.format(output_wksp.readY(0)[:5]))
print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5]))
Ouptut: Ouptut:
...@@ -74,6 +84,7 @@ Ouptut: ...@@ -74,6 +84,7 @@ Ouptut:
Input workspace: [ 38. 38. 38. 38. 38.] Input workspace: [ 38. 38. 38. 38. 38.]
Correction workspace: [ 10.26463773 9.81128219 9.39826191 9.02042771 8.67347109] 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: [ 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.** **Example - Basics of running CalculateEfficiencyCorrection with Density and ChemicalFormula.**
...@@ -86,17 +97,22 @@ Ouptut: ...@@ -86,17 +97,22 @@ Ouptut:
NumEvents=10000, NumBanks=1, BankPixelWidth=1) NumEvents=10000, NumBanks=1, BankPixelWidth=1)
# Calculate the efficiency correction # Calculate the efficiency correction
correction_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp,
Density=6.11, Density=6.11,
ChemicalFormula="V") 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 # Apply the efficiency correction to the measured spectrum
input_wksp = ConvertToPointData(InputWorkspace=input_wksp) 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('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: {}'.format(output_wksp.readY(0)[:5]))
print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5]))
Ouptut: Ouptut:
...@@ -105,6 +121,7 @@ Ouptut: ...@@ -105,6 +121,7 @@ Ouptut:
Input workspace: [ 38. 38. 38. 38. 38.] Input workspace: [ 38. 38. 38. 38. 38.]
Correction workspace: [ 24.40910309 23.29738394 22.28449939 21.35783225 20.50682528] 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: [ 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.** **Example - Basics of running CalculateEfficiencyCorrection with MeasuredEfficiency and ChemicalFormula.**
...@@ -117,17 +134,24 @@ Ouptut: ...@@ -117,17 +134,24 @@ Ouptut:
NumEvents=10000, NumBanks=1, BankPixelWidth=1) NumEvents=10000, NumBanks=1, BankPixelWidth=1)
# Calculate the efficiency correction # Calculate the efficiency correction
correction_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp, corr_wksp = CalculateEfficiencyCorrection(InputWorkspace=input_wksp,
MeasuredEfficiency=1e-2, MeasuredEfficiency=1e-2,
ChemicalFormula="(He3)") 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 # Apply the efficiency correction to the measured spectrum
input_wksp = ConvertToPointData(InputWorkspace=input_wksp) 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('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: {}'.format(output_wksp.readY(0)[:5]))
print('Output workspace using WavelengthRange: {}'.format(output_wksp_with_wave_range.readY(0)[:5]))
Ouptut: Ouptut:
...@@ -137,6 +161,83 @@ Ouptut: ...@@ -137,6 +161,83 @@ Ouptut:
Correction workspace: [ 873.27762699 832.68332786 795.69741128 761.85923269 730.78335476] Correction workspace: [ 873.27762699 832.68332786 795.69741128 761.85923269 730.78335476]
Output workspace: [ 33184.54982567 31641.9664586 30236.50162877 28950.65084207 Output workspace: [ 33184.54982567 31641.9664586 30236.50162877 28950.65084207
27769.76748099] 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.** **Example - Running CalculateEfficiencyCorrection for incident spectrum.**
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment