Skip to content

Commit

Permalink
Refactor more Python algorithms
Browse files Browse the repository at this point in the history
Refs #10706
  • Loading branch information
DanNixon committed Feb 17, 2015
1 parent 66e1452 commit b030a92
Show file tree
Hide file tree
Showing 6 changed files with 117 additions and 82 deletions.
Expand Up @@ -12,14 +12,17 @@ def category(self):


def summary(self):
return 'Calculates the scattering & transmission for a given sample material and size over a given wavelength range.'
return "Calculates the scattering & transmission for a given sample \
material and size over a given wavelength range."


def PyInit(self):
self.declareProperty(name='WavelengthRange', defaultValue='', validator=StringMandatoryValidator(),
self.declareProperty(name='WavelengthRange', defaultValue='',
validator=StringMandatoryValidator(),
doc='Wavelength range to calculate transmission for.')

self.declareProperty(name='ChemicalFormula', defaultValue='', validator=StringMandatoryValidator(),
self.declareProperty(name='ChemicalFormula', defaultValue='',
validator=StringMandatoryValidator(),
doc='Sample chemical formula')

self.declareProperty(name='NumberDensity', defaultValue=0.1,
Expand All @@ -29,18 +32,19 @@ def PyInit(self):
doc='Sample thickness (cm). Default=0.1')

self.declareProperty(MatrixWorkspaceProperty('OutputWorkspace', '', Direction.Output),
doc='Outputs the sample transmission over the wavelength range as a function of wavelength.')
doc='Outputs the sample transmission over the wavelength range '
'as a function of wavelength.')


def validateInputs(self):
issues = dict()

density = self.getProperty('NumberDensity').value
if(density < 0.0):
if density < 0.0:
issues['NumberDensity'] = 'NumberDensity must be positive'

thickness = self.getProperty('Thickness').value
if(thickness < 0.0):
if thickness < 0.0:
issues['Thickness'] = 'Thickness must be positive'

return issues
Expand All @@ -51,12 +55,13 @@ def PyExec(self):

# Create the workspace and set the sample material
CreateWorkspace(OutputWorkspace=self._output_ws, NSpec=2, DataX=[0, 0], DataY=[0, 0])
Rebin(InputWorkspace=self._output_ws, OutputWorkspace=self._output_ws, Params=self._bin_params)
Rebin(InputWorkspace=self._output_ws, OutputWorkspace=self._output_ws,
Params=self._bin_params)
SetSampleMaterial(InputWorkspace=self._output_ws, ChemicalFormula=self._chamical_formula)
ConvertToPointData(InputWorkspace=self._output_ws, OutputWorkspace=self._output_ws)

ws = mtd[self._output_ws]
wavelengths = ws.readX(0)
workspace = mtd[self._output_ws]
wavelengths = workspace.readX(0)
transmission_data = np.zeros(len(wavelengths))
scattering_data = np.zeros(len(wavelengths))

Expand All @@ -66,8 +71,8 @@ def PyExec(self):
transmission_data[idx] = transmission
scattering_data[idx] = scattering

ws.setY(0, transmission_data)
ws.setY(1, scattering_data)
workspace.setY(0, transmission_data)
workspace.setY(1, scattering_data)

self.setProperty('OutputWorkspace', self._output_ws)

Expand Down Expand Up @@ -103,6 +108,5 @@ def _calculate_at_wavelength(self, wavelength):
return transmission, scattering



# Register algorithm with Mantid
AlgorithmFactory.subscribe(CalculateSampleTransmission)
Expand Up @@ -5,11 +5,13 @@
class CreateEmptyTableWorkspace(PythonAlgorithm):

def summary(self):
return "Creates an empty TableWorkspace which can be populated with various types of information."
return "Creates an empty TableWorkspace which can be populated with various "\
"types of information."

def PyInit(self):
# Declare properties
self.declareProperty(ITableWorkspaceProperty("OutputWorkspace", "", Direction.Output), "The name of the table workspace that will be created.")
self.declareProperty(ITableWorkspaceProperty("OutputWorkspace", "", Direction.Output),
"The name of the table workspace that will be created.")

def PyExec(self):
tableWS = WorkspaceFactory.createTable()
Expand Down
Expand Up @@ -34,7 +34,7 @@
class LoadVesuvio(PythonAlgorithm):

def summary(self):
return "Loads raw data produced by the Vesuvio instrument at ISIS."
return "Loads raw data produced by the Vesuvio instrument at ISIS."

def PyInit(self):
self.declareProperty(RUN_PROP, "", StringMandatoryValidator(),
Expand All @@ -53,12 +53,13 @@ def PyInit(self):
self.declareProperty(FileProperty(INST_PAR_PROP,"",action=FileAction.OptionalLoad,
extensions=["dat"]),
doc="An optional IP file. If provided the values are used to correct "
"the default instrument values and attach the t0 values to each detector")
"the default instrument values and attach the t0 values to each "
"detector")

self.declareProperty(SUM_PROP, False,
doc="If true then the final output is a single spectrum containing the sum "
"of all of the requested spectra. All detector angles/parameters are "
"averaged over the individual inputs")
doc="If true then the final output is a single spectrum containing "
"the sum of all of the requested spectra. All detector angles/"
"parameters are averaged over the individual inputs")

self.declareProperty(WorkspaceProperty(WKSP_PROP, "", Direction.Output),
doc="The name of the output workspace.")
Expand Down Expand Up @@ -114,13 +115,13 @@ def _exec_single_foil_state_mode(self):
"""
runs = self._get_runs()
if len(runs) > 1:
raise RuntimeError("Single soil state mode does not currently support summing multiple files")
raise RuntimeError("Single soil state mode does not currently support summing "
"multiple files")

isis = config.getFacility("ISIS")
inst_prefix = isis.instrument("VESUVIO").shortName()

try:
run = int(runs[0])
run_str = inst_prefix + runs[0]
except ValueError:
run_str = runs[0]
Expand All @@ -138,7 +139,7 @@ def _exec_single_foil_state_mode(self):
foil_map = SpectraToFoilPeriodMap(self._nperiods)
for ws_index, spectrum_no in enumerate(all_spectra):
self._set_spectra_type(spectrum_no)
foil_out_periods, foil_thin_periods, foil_thick_periods = self._get_foil_periods()
foil_out_periods, foil_thin_periods, _ = self._get_foil_periods()

if self._diff_opt == "FoilOut":
raw_grp_indices = foil_map.get_indices(spectrum_no, foil_out_periods)
Expand Down Expand Up @@ -252,13 +253,15 @@ def _setup_raw(self, spectra):
# Cache delta_t values
raw_t = first_ws.readX(0)
delay = raw_t[2] - raw_t[1]
raw_t = raw_t - delay # The original EVS loader, raw.for/rawb.for, does this. Done here to match results
# The original EVS loader, raw.for/rawb.for, does this. Done here to match results
raw_t = raw_t - delay
self.pt_times = raw_t[1:]
self.delta_t = (raw_t[1:] - raw_t[:-1])

mon_raw_t = self._raw_monitors[0].readX(0)
delay = mon_raw_t[2] - mon_raw_t[1]
mon_raw_t = mon_raw_t - delay # The original EVS loader, raw.for/rawb.for, does this. Done here to match results
# The original EVS loader, raw.for/rawb.for, does this. Done here to match results
mon_raw_t = mon_raw_t - delay
self.mon_pt_times = mon_raw_t[1:]
self.delta_tmon = (mon_raw_t[1:] - mon_raw_t[:-1])

Expand Down Expand Up @@ -310,7 +313,8 @@ def _get_runs(self):
# Load is not doing the right thing when summing. The numbers don't look correct
if "-" in run_str:
lower,upper = run_str.split("-")
runs = range(int(lower),int(upper)+1) #range goes lower to up-1 but we want to include the last number
# Range goes lower to up-1 but we want to include the last number
runs = range(int(lower),int(upper)+1)

elif "," in run_str:
runs = run_str.split(",")
Expand All @@ -322,9 +326,9 @@ def _get_runs(self):

def _set_spectra_type(self, spectrum_no):
"""
Set whether this spectrum no is forward/backward scattering
and set the normalization range appropriately
@param spectrum_no The current spectrum no
Set whether this spectrum no is forward/backward scattering
and set the normalization range appropriately
@param spectrum_no The current spectrum no
"""
if spectrum_no >= self._backward_spectra_list[0] and spectrum_no <= self._backward_spectra_list[-1]:
self._spectra_type=BACKWARD
Expand Down Expand Up @@ -385,7 +389,8 @@ def _create_foil_workspaces(self):
nhists = first_ws.getNumberHistograms()
data_kwargs = {'NVectors':nhists,'XLength':ndata_bins,'YLength':ndata_bins}

self.foil_out = WorkspaceFactory.create(first_ws, **data_kwargs) # This will be used as the result workspace
# This will be used as the result workspace
self.foil_out = WorkspaceFactory.create(first_ws, **data_kwargs)
self.foil_out.setDistribution(True)
self.foil_thin = WorkspaceFactory.create(first_ws, **data_kwargs)

Expand Down Expand Up @@ -417,20 +422,21 @@ def _sum_foil_periods(self):
foil_out_periods, foil_thin_periods, foil_thick_periods = self._get_foil_periods()

if self._nperiods == 6 and self._spectra_type == FORWARD:
mon_out_periods = (5,6)
mon_thin_periods = (3,4)
mon_thick_periods = foil_thick_periods
mon_out_periods = (5,6)
mon_thin_periods = (3,4)
mon_thick_periods = foil_thick_periods
else:
# None indicates same as standard foil
mon_out_periods, mon_thin_periods, mon_thick_periods = (None,None,None)
#

# Foil out
self._sum_foils(self.foil_out, self.mon_out, IOUT, foil_out_periods, mon_out_periods)
# Thin foil
self._sum_foils(self.foil_thin, self.mon_thin, ITHIN, foil_thin_periods, mon_thin_periods)
# Thick foil
if foil_thick_periods is not None:
self._sum_foils(self.foil_thick, self.mon_thick, ITHICK, foil_thick_periods, mon_thick_periods)
self._sum_foils(self.foil_thick, self.mon_thick, ITHICK,
foil_thick_periods, mon_thick_periods)

#----------------------------------------------------------------------------------------
def _get_foil_periods(self):
Expand Down Expand Up @@ -467,13 +473,14 @@ def _get_foil_periods(self):
#----------------------------------------------------------------------------------------
def _sum_foils(self, foil_ws, mon_ws, sum_index, foil_periods, mon_periods=None):
"""
Sums the counts from the given foil periods in the raw data group
@param foil_ws :: The workspace that will receive the summed counts
@param mon_ws :: The monitor workspace that will receive the summed monitor counts
@param sum_index :: An index into the sum3 array where the integrated counts will be accumulated
@param foil_periods :: The period numbers that contribute to this sum
@param mon_periods :: The period numbers of the monitors that contribute to this monitor sum
(if None then uses the foil_periods)
Sums the counts from the given foil periods in the raw data group
@param foil_ws :: The workspace that will receive the summed counts
@param mon_ws :: The monitor workspace that will receive the summed monitor counts
@param sum_index :: An index into the sum3 array where the integrated counts will be
accumulated
@param foil_periods :: The period numbers that contribute to this sum
@param mon_periods :: The period numbers of the monitors that contribute to this monitor sum
(if None then uses the foil_periods)
"""
raw_grp_indices = self.foil_map.get_indices(self._spectrum_no, foil_periods)
wsindex = self._ws_index
Expand Down Expand Up @@ -511,7 +518,9 @@ def _normalise_by_monitor(self):
wsindex = self._ws_index
# inner function to apply normalization
def monitor_normalization(foil_ws, mon_ws):
"""Applies monitor normalization to the given foil spectrum from the given monitor spectrum
"""
Applies monitor normalization to the given foil spectrum from the given
monitor spectrum.
"""
mon_values = mon_ws.readY(wsindex)
mon_values_sum = np.sum(mon_values[indices_in_range])
Expand Down Expand Up @@ -542,7 +551,8 @@ def normalise_to_out(foil_ws, foil_type):
values = foil_ws.dataY(wsindex)
sum_values = np.sum(values[range_indices])
if sum_values == 0.0:
self.getLogger().warning("No counts in %s foil spectrum %d." % (foil_type,self._spectrum_no))
self.getLogger().warning("No counts in %s foil spectrum %d." % (
foil_type,self._spectrum_no))
sum_values = 1.0
norm_factor = (sum_out/sum_values)
values *= norm_factor
Expand Down Expand Up @@ -611,7 +621,8 @@ def _calculate_double_difference(self, ws_index):
eout = self.foil_out.dataE(ws_index)
ethin = self.foil_thin.readE(ws_index)
ethick = self.foil_thick.readE(ws_index)
np.sqrt((one_min_beta*eout)**2 + ethin**2 + (self._beta**2)*ethick**2, eout) # The second argument makes it happen in place
# The second argument makes it happen in place
np.sqrt((one_min_beta*eout)**2 + ethin**2 + (self._beta**2)*ethick**2, eout)

#----------------------------------------------------------------------------------------
def _calculate_thick_difference(self, ws_index):
Expand Down Expand Up @@ -693,7 +704,8 @@ def _get_header_format(self, ip_filename):
try:
return IP_HEADERS[len(columns)]
except KeyError:
raise ValueError("Unknown format for IP file. Currently support 5/6 column variants. ncols=%d" % (len(columns)))
raise ValueError("Unknown format for IP file. Currently support 5/6 column "
"variants. ncols=%d" % (len(columns)))

#----------------------------------------------------------------------------------------
def _store_results(self):
Expand Down Expand Up @@ -794,11 +806,13 @@ def get_foil_periods(self, spectrum_no, state):
return foil_periods

def get_indices(self, spectrum_no, foil_state_numbers):
"""Returns a tuple of indices that can be used to access the Workspace within
"""
Returns a tuple of indices that can be used to access the Workspace within
a WorkspaceGroup that corresponds to the foil state numbers given
@param spectrum_no :: A spectrum number (1->nspectra)
@param foil_state_no :: A number between 1 & 9(inclusive) that defines which foil state is required
@returns A tuple of indices in a WorkspaceGroup that gives the associated Workspace
@param spectrum_no :: A spectrum number (1->nspectra)
@param foil_state_no :: A number between 1 & 9(inclusive) that defines which foil
state is required
@returns A tuple of indices in a WorkspaceGroup that gives the associated Workspace
"""
indices = []
for state in foil_state_numbers:
Expand Down Expand Up @@ -834,11 +848,13 @@ def get_index(self, spectrum_no, foil_state_no):

def _validate_foil_number(self, foil_number):
if foil_number < 1 or foil_number > 9:
raise ValueError("Invalid foil state given, expected a number between 1 and 9. number=%d" % foil_number)
raise ValueError("Invalid foil state given, expected a number between "
"1 and 9. number=%d" % foil_number)

def _validate_spectrum_number(self, spectrum_no):
if spectrum_no < 1 or spectrum_no > 198:
raise ValueError("Invalid spectrum given, expected a number between 3 and 198. spectrum=%d" % spectrum_no)
raise ValueError("Invalid spectrum given, expected a number between 3 "
"and 198. spectrum=%d" % spectrum_no)

#########################################################################################

Expand Down
@@ -1,9 +1,8 @@
from mantid.api import PythonAlgorithm, AlgorithmFactory,WorkspaceProperty,PropertyMode
from mantid.kernel import Direction,IntArrayProperty, FloatArrayProperty
from mantid.api import PythonAlgorithm, AlgorithmFactory
from mantid.kernel import Direction, IntArrayProperty, FloatArrayProperty
import mantid,math,numpy



class SortDetectors(PythonAlgorithm):
""" Sort detectors by distance
"""
Expand All @@ -25,30 +24,32 @@ def summary(self):
def PyInit(self):
""" Declare properties
"""
self.declareProperty(mantid.api.WorkspaceProperty("Workspace","",direction=mantid.kernel.Direction.Input, validator=mantid.api.InstrumentValidator()), "Input workspace")
self.declareProperty(mantid.api.WorkspaceProperty("Workspace", "",
direction=mantid.kernel.Direction.Input,
validator=mantid.api.InstrumentValidator()),
"Input workspace")

self.declareProperty(IntArrayProperty("UpstreamSpectra",Direction.Output))
self.declareProperty(FloatArrayProperty("UpstreamDetectorDistances",Direction.Output))
self.declareProperty(IntArrayProperty("DownstreamSpectra",Direction.Output))
self.declareProperty(FloatArrayProperty("DownstreamDetectorDistances",Direction.Output))
return
self.declareProperty(IntArrayProperty("UpstreamSpectra", Direction.Output))
self.declareProperty(FloatArrayProperty("UpstreamDetectorDistances", Direction.Output))
self.declareProperty(IntArrayProperty("DownstreamSpectra", Direction.Output))
self.declareProperty(FloatArrayProperty("DownstreamDetectorDistances", Direction.Output))

def PyExec(self):
""" Main execution body
"""
w = self.getProperty("Workspace").value
samplePos=w.getInstrument().getSample().getPos()
moderatorPos=w.getInstrument().getSource().getPos()
incident=samplePos-moderatorPos
workspace = self.getProperty("Workspace").value
samplePos = workspace.getInstrument().getSample().getPos()
moderatorPos = workspace.getInstrument().getSource().getPos()
incident = samplePos - moderatorPos

upstream=[]
upinds=[]
updist=[]
downstream=[]
downinds=[]
downdist=[]
for i in range(w.getNumberHistograms()):
detPos=w.getDetector(i).getPos()
for i in range(workspace.getNumberHistograms()):
detPos=workspace.getDetector(i).getPos()
scattered=detPos-samplePos
if abs(scattered.angle(incident))>0.999*math.pi:
upstream.append((i,scattered.norm()))
Expand Down
Expand Up @@ -111,8 +111,9 @@ def _calculate_resolution(self, workspace, output_ws_name):
fit.initialize()
fit.setChild(True)
mantid.simpleapi._set_properties(fit, function, InputWorkspace=workspace, MaxIterations=0,
CreateOutput=True, Output=fit_naming_stem, WorkspaceIndex=self._spectrum_index,
OutputCompositeMembers=True)
CreateOutput=True, Output=fit_naming_stem,
WorkspaceIndex=self._spectrum_index,
OutputCompositeMembers=True)
fit.execute()
fit_ws = fit.getProperty('OutputWorkspace').value

Expand Down

0 comments on commit b030a92

Please sign in to comment.