Skip to content

Commit

Permalink
More refactoring of algorithm
Browse files Browse the repository at this point in the history
Refs #1118
  • Loading branch information
DanNixon committed Feb 27, 2015
1 parent 82acd04 commit a84ca45
Showing 1 changed file with 85 additions and 78 deletions.
Original file line number Diff line number Diff line change
@@ -1,33 +1,35 @@
from mantid.simpleapi import *
from mantid.api import PythonAlgorithm, AlgorithmFactory, MatrixWorkspaceProperty, WorkspaceGroupProperty
from mantid.api import PythonAlgorithm, AlgorithmFactory, PropertyMode, StringMandatoryValidator,\
MatrixWorkspaceProperty, WorkspaceGroupProperty
from mantid.kernel import StringListValidator, Direction, logger
import math, numpy as np


class FlatPaalmanPingsCorrection(PythonAlgorithm):
class FlatPlatePaalmanPingsCorrection(PythonAlgorithm):

_sample_ws_name = None
_sample_chemical_formula = None
_sample_number_density = None
_sample_thickness = None
_sample_angle = 0.0
_usecan = False
_use_can = False
_can_ws_name = None
_can_chemical_formula = None
_can_number_density = None
_can_thickness1 = None
_can_thickness2 = None
_can_front_thickness = None
_can_back_thickness = None
_can_scale = 1.0
_number_wavelengths = 10
_emode = None
_efixed = 0.0
_interpolate = 'Linear'
_output_ws_name = None
_angles = list()
_waves = list()
_elastic = 0.0


def category(self):
return "Workflow\\MIDAS;PythonAlgorithms"
return "Workflow\\MIDAS;PythonAlgorithms;CorrectionFunctions\\AbsorptionCorrections"


def summary(self):
Expand All @@ -40,21 +42,23 @@ def PyInit(self):
doc='Name for the input Sample workspace.')

self.declareProperty(name='SampleChemicalFormula', defaultValue='',
validator=StringMandatoryValidator(),
doc='Sample chemical formula')
self.declareProperty(name='SampleNumberDensity', defaultValue=0.0,
self.declareProperty(name='SampleNumberDensity', defaultValue=0.1,
doc='Sample number density')
self.declareProperty(name='SampleThickness', defaultValue=0.0,
doc='Sample thickness')
self.declareProperty(name='SampleAngle', defaultValue=0.0,
doc='Sample angle')

self.declareProperty(MatrixWorkspaceProperty('CanWorkspace', '',
direction=Direction.Input),
direction=Direction.Input,
optional=PropertyMode.Optional),
doc="Name for the input Can workspace.")

self.declareProperty(name='CanChemicalFormula', defaultValue='',
doc='Can chemical formula')
self.declareProperty(name='CanNumberDensity', defaultValue=0.0,
self.declareProperty(name='CanNumberDensity', defaultValue=0.1,
doc='Can number density')

self.declareProperty(name='CanFrontThickness', defaultValue=0.0,
Expand Down Expand Up @@ -87,63 +91,64 @@ def PyExec(self):
SampleNumberDensity=self._sample_number_density)

# If using a can, set sample material using chemical formula
if self._usecan:
if self._use_can:
SetSampleMaterial(InputWorkspace=self._can_ws_name, ChemicalFormula=self._can_chemical_formula,
SampleNumberDensity=self._can_number_density)

dataA1 = []
dataA2 = []
dataA3 = []
dataA4 = []
# Holders for the corrected data
data_ass = []
data_assc = []
data_acsc = []
data_acc = []

self._get_angles()
num_angles = len(self._angles)

for angle_idx in range(num_angles):
angle = self._angles[angle_idx]
(A1, A2, A3, A4) = self._flatAbs(angle)
(ass, assc, acsc, acc) = self._flat_abs(angle)

logger.information('Angle ' + str(angle_idx+1) + ' : ' + str(self._angles[angle_idx]) + ' successful')
logger.information('Angle %d: %f successful' % (angle_idx+1, self._angles[angle_idx]))

dataA1 = np.append(dataA1, A1)
dataA2 = np.append(dataA2, A2)
dataA3 = np.append(dataA3, A3)
dataA4 = np.append(dataA4, A4)
data_ass = np.append(data_ass, ass)
data_assc = np.append(data_assc, assc)
data_acsc = np.append(data_acsc, acsc)
data_acc = np.append(data_acc, acc)

sample_logs = {'sample_shape': 'flatplate', 'sample_filename': self._sample_ws_name,
'sample_thickness': self._sample_thickness, 'sample_angle': self._sample_angle}
dataX = self._waves * num_angles

# Create the output workspaces
ass_ws = self._output_ws_name + '_ass'
CreateWorkspace(OutputWorkspace=ass_ws, DataX=dataX, DataY=dataA1,
CreateWorkspace(OutputWorkspace=ass_ws, DataX=dataX, DataY=data_ass,
NSpec=num_angles, UnitX='Wavelength')
self._addSampleLogs(ass_ws, sample_logs)
self._add_sample_logs(ass_ws, sample_logs)

workspaces = [ass_ws]

if self._usecan:
if self._use_can:
AddSampleLog(Workspace=ass_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))

assc_ws = self._output_ws_name + '_assc'
workspaces.append(assc_ws)
CreateWorkspace(OutputWorkspace=assc_ws, DataX=dataX, DataY=dataA2,
CreateWorkspace(OutputWorkspace=assc_ws, DataX=dataX, DataY=data_assc,
NSpec=num_angles, UnitX='Wavelength')
self._addSampleLogs(assc_ws, sample_logs)
self._add_sample_logs(assc_ws, sample_logs)
AddSampleLog(Workspace=assc_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))

acsc_ws = self._output_ws_name + '_acsc'
workspaces.append(acsc_ws)
CreateWorkspace(OutputWorkspace=acsc_ws, DataX=dataX, DataY=dataA3,
CreateWorkspace(OutputWorkspace=acsc_ws, DataX=dataX, DataY=data_acsc,
NSpec=num_angles, UnitX='Wavelength')
self._addSampleLogs(acsc_ws, sample_logs)
self._add_sample_logs(acsc_ws, sample_logs)
AddSampleLog(Workspace=acsc_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))

acc_ws = self._output_ws_name + '_acc'
workspaces.append(acc_ws)
CreateWorkspace(OutputWorkspace=acc_ws, DataX=dataX, DataY=dataA4,
CreateWorkspace(OutputWorkspace=acc_ws, DataX=dataX, DataY=data_acc,
NSpec=num_angles, UnitX='Wavelength')
self._addSampleLogs(acc_ws, sample_logs)
self._add_sample_logs(acc_ws, sample_logs)
AddSampleLog(Workspace=acc_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))

GroupWorkspaces(InputWorkspaces=','.join(workspaces), OutputWorkspace=self._output_ws_name)
Expand All @@ -158,12 +163,12 @@ def _setup(self):
self._sample_angle = self.getProperty('SampleAngle').value

self._can_ws_name = self.getPropertyValue('CanWorkspace')
self._usecan = self._can_ws_name != ''
self._use_can = self._can_ws_name != ''

self._can_chemical_formula = self.getPropertyValue('CanChemicalFormula')
self._can_number_density = self.getProperty('CanNumberDensity').value
self._can_thickness1 = self.getProperty('CanFrontThickness').value
self._can_thickness2 = self.getProperty('CanBackThickness').value
self._can_front_thickness = self.getProperty('CanFrontThickness').value
self._can_back_thickness = self.getProperty('CanBackThickness').value
self._can_scale = self.getProperty('CanScaleFactor').value

self._number_wavelengths = self.getProperty('NumberWavelengths').value
Expand Down Expand Up @@ -191,20 +196,21 @@ def _wave_range(self):
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)
waves = []
for l in range(0,number_waves):
waves.append(wave_min+l*wave_bin)
self._waves = waves
wave_bin = (wave_max - wave_min) / (number_waves-1)

for idx in range(0, number_waves):
self._waves.append(wave_min + idx * wave_bin)

if self._emode == 'Elastic':
self._elastic = waves[int(number_waves / 2)]
if self._emode == 'Indirect':
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)
DeleteWorkspace(wave_range)


def _addSampleLogs(self, ws, sample_logs):
def _add_sample_logs(self, ws, sample_logs):
"""
Add a dictionary of logs to a workspace.
Expand All @@ -225,16 +231,15 @@ def _addSampleLogs(self, ws, sample_logs):
AddSampleLog(Workspace=ws, LogName=key, LogType=log_type, LogText=str(value))


def _flatAbs(self, angle):
def _flat_abs(self, angle):
"""
FlatAbs - calculate flat plate absorption factors
For more information See:
- MODES User Guide: http://www.isis.stfc.ac.uk/instruments/iris/data-analysis/modes-v3-user-guide-6962.pdf
- C J Carlile, Rutherford Laboratory report, RL-74-103 (1974)
@param angle - angle
"""

PICONV = math.pi / 180.0

canAngle = self._sample_angle * PICONV
Expand Down Expand Up @@ -263,36 +268,35 @@ def _flatAbs(self, angle):
sec1 = 1.0 / math.cos(canAngle)
sec2 = 1.0 / math.cos(tsec)

#list of wavelengths
# List of wavelengths
waves = np.array(self._waves)

#sample cross section
sampleXSection = (sam_material.totalScatterXSection() + sam_material.absorbXSection() * waves / 1.8) * self._sample_number_density
# Sample cross section
sample_x_section = (sam_material.totalScatterXSection() + sam_material.absorbXSection() * waves / 1.8) * self._sample_number_density

#vector version of fact
vecFact = np.vectorize(self._fact)
fs = vecFact(sampleXSection, self._sample_thickness, sec1, sec2)
fs = vecFact(sample_x_section, self._sample_thickness, sec1, sec2)

sampleSec1, sampleSec2 = self._calc_thickness_at_x_sect(sampleXSection, self._sample_thickness, [sec1, sec2])
sample_sect_1, sample_sect_2 = self._calc_thickness_at_x_sect(sample_x_section, self._sample_thickness, [sec1, sec2])

if sec2 < 0.0:
ass = fs / self._sample_thickness
else:
ass= np.exp(-sampleSec2) * fs / self._sample_thickness
ass= np.exp(-sample_sect_2) * fs / self._sample_thickness

if self._usecan:
if self._use_can:
can_sample = mtd[self._can_ws_name].sample()
can_material = can_sample.getMaterial()

#calculate can cross section
canXSection = (can_material.totalScatterXSection() + can_material.absorbXSection() * waves / 1.8) * self._can_number_density
assc, acsc, acc = self._calculate_can(ass, canXSection, self._can_thickness1, self._can_thickness2, sampleSec1, sampleSec2, [sec1, sec2])
# Calculate can cross section
can_x_section = (can_material.totalScatterXSection() + can_material.absorbXSection() * waves / 1.8) * self._can_number_density
assc, acsc, acc = self._calculate_can(ass, can_x_section, sample_sect_1, sample_sect_2, [sec1, sec2])

return ass, assc, acsc, acc


def _fact(self, xSection, thickness, sec1, sec2):
S = xSection * thickness * (sec1 - sec2)
def _fact(self, x_section, thickness, sec1, sec2):
S = x_section * thickness * (sec1 - sec2)
F = 1.0
if S == 0.0:
F = thickness
Expand All @@ -302,58 +306,61 @@ def _fact(self, xSection, thickness, sec1, sec2):
return F


def _calc_thickness_at_x_sect(self, xSection, thickness, sec):
def _calc_thickness_at_x_sect(self, x_section, thickness, sec):
sec1, sec2 = sec

thickSec1 = xSection * thickness * sec1
thickSec2 = xSection * thickness * sec2
thick_sec_1 = x_section * thickness * sec1
thick_sec_2 = x_section * thickness * sec2

return thick_sec_1, thick_sec_2

return thickSec1, thickSec2

def _calculate_can(self, ass, can_x_section, sample_sect_1, sample_sect_2, sec):
"""
Calculates the A_s,sc, A_c,sc and A_c,c data.
"""

def _calculate_can(self, ass, canXSection, canThickness1, canThickness2, sampleSec1, sampleSec2, sec):
assc = np.ones(ass.size)
acsc = np.ones(ass.size)
acc = np.ones(ass.size)

sec1, sec2 = sec

#vector version of fact
vecFact = np.vectorize(self._fact)
f1 = vecFact(canXSection, canThickness1, sec1, sec2)
f2 = vecFact(canXSection, canThickness2, sec1, sec2)
f1 = vecFact(can_x_section, self._can_front_thickness, sec1, sec2)
f2 = vecFact(can_x_section, self._can_back_thickness, sec1, sec2)

canThick1Sec1, canThick1Sec2 = self._calc_thickness_at_x_sect(canXSection, canThickness1, sec)
_, canThick2Sec2 = self._calc_thickness_at_x_sect(canXSection, canThickness2, sec)
can_thick_1_sect_1, can_thick_1_sect_2 = self._calc_thickness_at_x_sect(can_x_section, self._can_front_thickness, sec)
_, can_thick_2_sect_2 = self._calc_thickness_at_x_sect(can_x_section, self._can_back_thickness, sec)

if sec2 < 0.0:
val = np.exp(-(canThick1Sec1-canThick1Sec2))
val = np.exp(-(can_thick_1_sect_1 - can_thick_1_sect_2))
assc = ass * val

acc1 = f1
acc2 = f2 * val

acsc1 = acc1
acsc2 = acc2 * np.exp(-(sampleSec1 - sampleSec2))
acsc2 = acc2 * np.exp(-(sample_sect_1 - sample_sect_2))

else:
val = np.exp(-(canThick1Sec1 + canThick2Sec2))
val = np.exp(-(can_thick_1_sect_1 + can_thick_2_sect_2))
assc = ass * val

acc1 = f1 * np.exp(-(canThick1Sec2 + canThick2Sec2))
acc1 = f1 * np.exp(-(can_thick_1_sect_2 + can_thick_2_sect_2))
acc2 = f2 * val

acsc1 = acc1 * np.exp(-sampleSec2)
acsc2 = acc2 * np.exp(-sampleSec1)
acsc1 = acc1 * np.exp(-sample_sect_2)
acsc2 = acc2 * np.exp(-sample_sect_1)

canThickness = canThickness1 + canThickness2
can_thickness = self._can_front_thickness + self._can_back_thickness

if canThickness > 0.0:
acc = (acc1 + acc2) / canThickness
acsc = (acsc1 + acsc2) / canThickness
if can_thickness > 0.0:
acc = (acc1 + acc2) / can_thickness
acsc = (acsc1 + acsc2) / can_thickness

return assc, acsc, acc


# Register algorithm with Mantid
AlgorithmFactory.subscribe(FlatPaalmanPingsCorrection)
AlgorithmFactory.subscribe(FlatPlatePaalmanPingsCorrection)

0 comments on commit a84ca45

Please sign in to comment.