diff --git a/Framework/PythonInterface/plugins/algorithms/BASISReduction.py b/Framework/PythonInterface/plugins/algorithms/BASISReduction.py index e94b1c7702115480997e331edc79f298e288e42c..35b7c4ee2b4a1d44666c65ae3f44ab6011c057d2 100644 --- a/Framework/PythonInterface/plugins/algorithms/BASISReduction.py +++ b/Framework/PythonInterface/plugins/algorithms/BASISReduction.py @@ -4,8 +4,10 @@ from __future__ import (absolute_import, division, print_function) import numpy as np import mantid.simpleapi as sapi -from mantid.api import mtd, PythonAlgorithm, AlgorithmFactory, FileProperty, FileAction -from mantid.kernel import IntArrayProperty, StringListValidator, FloatArrayProperty, EnabledWhenProperty,\ +from mantid.api import mtd, PythonAlgorithm, AlgorithmFactory,\ + FileProperty, FileAction +from mantid.kernel import IntArrayProperty, StringListValidator,\ + FloatArrayProperty, EnabledWhenProperty,\ FloatArrayLengthValidator, Direction, PropertyCriterion from mantid import config from os.path import join as pjoin @@ -106,7 +108,8 @@ class BASISReduction(PythonAlgorithm): self.declareProperty("GroupDetectors", "None", StringListValidator(grouping_type), "Switch for grouping detectors") - self.declareProperty("NormalizeToFirst", False, "Normalize spectra to intensity of spectrum with lowest Q?") + self.declareProperty("NormalizeToFirst", False, "Normalize spectra \ + to intensity of spectrum with lowest Q?") # Properties affected by the reflection selected titleReflection = "Reflection Selector" @@ -114,7 +117,8 @@ class BASISReduction(PythonAlgorithm): default_reflection = REFLECTIONS_DICT["silicon111"] self.declareProperty("ReflectionType", default_reflection["name"], StringListValidator(available_reflections), - "Analyzer. Documentation lists typical associated property values.") + "Analyzer. Documentation lists typical \ + associated property values.") self.setPropertyGroup("ReflectionType", titleReflection) self.declareProperty(FloatArrayProperty("EnergyBins", default_reflection["energy_bins"], @@ -127,7 +131,8 @@ class BASISReduction(PythonAlgorithm): "Momentum transfer binning scheme") self.setPropertyGroup("MomentumTransferBins", titleReflection) self.declareProperty(FileProperty(name="MaskFile", defaultValue='', - action=FileAction.OptionalLoad, extensions=['.xml']), + action=FileAction.OptionalLoad, + extensions=['.xml']), "See documentation for latest mask files.") self.setPropertyGroup("MaskFile", titleReflection) @@ -150,8 +155,10 @@ class BASISReduction(PythonAlgorithm): self.setPropertySettings("NormRunNumbers", ifDivideByVanadium) self.setPropertyGroup("NormRunNumbers", titleDivideByVanadium) arrVal = FloatArrayLengthValidator(2) - self.declareProperty(FloatArrayProperty("NormWavelengthRange", DEFAULT_RANGE, - arrVal, direction=Direction.Input), + self.declareProperty(FloatArrayProperty("NormWavelengthRange", + DEFAULT_RANGE, + arrVal, + direction=Direction.Input), "Wavelength range for normalization") self.setPropertySettings("NormWavelengthRange", ifDivideByVanadium) self.setPropertyGroup("NormWavelengthRange", titleDivideByVanadium) @@ -164,17 +171,20 @@ class BASISReduction(PythonAlgorithm): PropertyCriterion.IsNotDefault) self.setPropertyGroup('SaveNXSPE', title_nxspe) self.declareProperty('PsiAngleLog', 'SE50Rot', direction=Direction.Input, - doc='name of entry in the logs storing the psi angle') + doc='name of entry in the logs storing the psi \ + angle') self.setPropertySettings('PsiAngleLog', nxspe_enabled) self.setPropertyGroup('PsiAngleLog', title_nxspe) self.declareProperty('PsiOffset', 0.0, direction=Direction.Input, - doc='add this quantity to the psi angle stored in the log') + doc='add this quantity to the psi angle stored \ + in the log') self.setPropertySettings('PsiOffset', nxspe_enabled) self.setPropertyGroup('PsiOffset', title_nxspe) # Aditional output properties titleAddionalOutput = "Additional Output" - self.declareProperty("OutputSusceptibility", False, direction=Direction.Input, + self.declareProperty("OutputSusceptibility", False, + direction=Direction.Input, doc="Output dynamic susceptibility (Xqw)") self.setPropertyGroup("OutputSusceptibility", titleAddionalOutput) @@ -182,17 +192,19 @@ class BASISReduction(PythonAlgorithm): def PyExec(self): config['default.facility'] = "SNS" config['default.instrument'] = self._long_inst - self._reflection = REFLECTIONS_DICT[self.getProperty("ReflectionType").value] + self._reflection =\ + REFLECTIONS_DICT[self.getProperty("ReflectionType").value] self._doIndiv = self.getProperty("DoIndividual").value - self._etBins = 1.E-03 * self.getProperty("EnergyBins").value # micro-eV to mili-eV + # micro-eV to mili-eV + self._etBins = 1.E-03 * self.getProperty("EnergyBins").value self._qBins = self.getProperty("MomentumTransferBins").value - self._qBins[0] -= self._qBins[1]/2.0 # self._qBins[0] is leftmost bin boundary - self._qBins[2] += self._qBins[1]/2.0 # self._qBins[2] is rightmost bin boundary + self._qBins[0] -= self._qBins[1]/2.0 # leftmost bin boundary + self._qBins[2] += self._qBins[1]/2.0 # rightmost bin boundary self._noMonNorm = self.getProperty("NoMonitorNorm").value self._maskFile = self.getProperty("MaskFile").value maskfile = self.getProperty("MaskFile").value - self._maskFile = maskfile if maskfile else pjoin(DEFAULT_MASK_GROUP_DIR, - self._reflection["mask_file"]) + self._maskFile = maskfile if maskfile else\ + pjoin(DEFAULT_MASK_GROUP_DIR, self._reflection["mask_file"]) self._groupDetOpt = self.getProperty("GroupDetectors").value self._normalizeToFirst = self.getProperty("NormalizeToFirst").value self._doNorm = self.getProperty("DivideByVanadium").value @@ -213,7 +225,9 @@ class BASISReduction(PythonAlgorithm): config.appendDataSearchDir(DEFAULT_MASK_GROUP_DIR) self._maskFile = self._reflection["mask_file"] - sapi.LoadMask(Instrument='BASIS', OutputWorkspace='BASIS_MASK', InputFile=self._maskFile) + sapi.LoadMask(Instrument='BASIS', + OutputWorkspace='BASIS_MASK', + InputFile=self._maskFile) # Work around length issue _dMask = sapi.ExtractMask('BASIS_MASK') @@ -229,49 +243,63 @@ class BASISReduction(PythonAlgorithm): if ";" in norm_runs: raise SyntaxError("Normalization does not support run groups") self._normalizationType = self.getProperty("NormalizationType").value - self.log().information("Divide by Vanadium with normalization" + self._normalizationType) + self.log().information("Divide by Vanadium with normalization" + + self._normalizationType) - # The following steps are common to all types of Vanadium normalization + # Following steps common to all types of Vanadium normalization # norm_runs encompasses a single set, thus _getRuns returns # a list of only one item norm_set = self._getRuns(norm_runs, doIndiv=False)[0] normWs = self._sum_and_calibrate(norm_set, extra_extension="_norm") + normRange = self.getProperty("NormWavelengthRange").value + bin_width = normRange[1] - normRange[0] # This rebin integrates counts onto a histogram of a single bin - if self._normalizationType == "by detectorID": - normRange = self.getProperty("NormWavelengthRange").value - self._normRange = [normRange[0], normRange[1]-normRange[0], normRange[1]] - sapi.Rebin(InputWorkspace=normWs, OutputWorkspace=normWs, Params=self._normRange) - + if self._normalizationType == "by detector ID": + self._normRange = [normRange[0], bin_width, normRange[1]] + sapi.Rebin(InputWorkspace=normWs, + OutputWorkspace=normWs, + Params=self._normRange) + self._normWs = normWs # FindDetectorsOutsideLimits to be substituted by MedianDetectorTest - sapi.FindDetectorsOutsideLimits(InputWorkspace=normWs, OutputWorkspace="BASIS_NORM_MASK") - + sapi.FindDetectorsOutsideLimits(InputWorkspace=normWs, + LowThreshold=1.0*bin_width, + # no count events outside ranges + RangeLower=normRange[0], + RangeUpper=normRange[1], + OutputWorkspace='BASIS_NORM_MASK') # additional reduction steps when normalizing by Q slice if self._normalizationType == "by Q slice": - self._normWs = self._group_and_SofQW(normWs, self._etBins, isSample=False) - if not self._debugMode: - sapi.DeleteWorkspace(normWs) # Delete vanadium events file + self._normWs = self._group_and_SofQW(normWs, self._etBins, + isSample=False) ########################## ## Process the sample ## ########################## - self._run_list = self._getRuns(self.getProperty("RunNumbers").value, doIndiv=self._doIndiv) + self._run_list = self._getRuns(self.getProperty("RunNumbers").value, + doIndiv=self._doIndiv) for run_set in self._run_list: self._samWs = self._sum_and_calibrate(run_set) self._samWsRun = str(run_set[0]) # Divide by Vanadium detector ID, if pertinent if self._normalizationType == "by detector ID": # Mask detectors with insufficient Vanadium signal before dividing - sapi.MaskDetectors(Workspace=self._samWs, MaskedWorkspace='BASIS_NORM_MASK') - sapi.Divide(LHSWorkspace=self._samWs, RHSWorkspace=self._normWs, OutputWorkspace=self._samWs) + sapi.MaskDetectors(Workspace=self._samWs, + MaskedWorkspace='BASIS_NORM_MASK') + sapi.Divide(LHSWorkspace=self._samWs, + RHSWorkspace=self._normWs, + OutputWorkspace=self._samWs) # additional reduction steps - self._samSqwWs = self._group_and_SofQW(self._samWs, self._etBins, isSample=True) + self._samSqwWs = self._group_and_SofQW(self._samWs, self._etBins, + isSample=True) if not self._debugMode: sapi.DeleteWorkspace(self._samWs) # delete events file # Divide by Vanadium Q slice, if pertinent if self._normalizationType == "by Q slice": - sapi.Divide(LHSWorkspace=self._samSqwWs, RHSWorkspace=self._normWs, OutputWorkspace=self._samSqwWs) + sapi.Divide(LHSWorkspace=self._samSqwWs, + RHSWorkspace=self._normWs, + OutputWorkspace=self._samSqwWs) # Clear mask from reduced file. Needed for binary operations # involving this S(Q,w) sapi.ClearMaskFlag(Workspace=self._samSqwWs) @@ -280,48 +308,65 @@ class BASISReduction(PythonAlgorithm): self._ScaleY(self._samSqwWs) # Transform the vertical axis (Q) to point data + # Q-values are in X-axis now sapi.Transpose(InputWorkspace=self._samSqwWs, - OutputWorkspace=self._samSqwWs) # Q-values are in X-axis now + OutputWorkspace=self._samSqwWs) + # from histo to point sapi.ConvertToPointData(InputWorkspace=self._samSqwWs, - OutputWorkspace=self._samSqwWs) # from histo to point + OutputWorkspace=self._samSqwWs) + # Q-values back to vertical axis sapi.Transpose(InputWorkspace=self._samSqwWs, - OutputWorkspace=self._samSqwWs) # Q-values back to vertical axis + OutputWorkspace=self._samSqwWs) # Output Dave and Nexus files extension = "_divided.dat" if self._doNorm else ".dat" - dave_grp_filename = self._makeRunName(self._samWsRun, False) + extension - sapi.SaveDaveGrp(Filename=dave_grp_filename, InputWorkspace=self._samSqwWs, ToMicroEV=True) + dave_grp_filename = self._makeRunName(self._samWsRun, False) +\ + extension + sapi.SaveDaveGrp(Filename=dave_grp_filename, + InputWorkspace=self._samSqwWs, + ToMicroEV=True) extension = "_divided_sqw.nxs" if self._doNorm else "_sqw.nxs" - processed_filename = self._makeRunName(self._samWsRun, False) + extension - sapi.SaveNexus(Filename=processed_filename, InputWorkspace=self._samSqwWs) + processed_filename = self._makeRunName(self._samWsRun, False) +\ + extension + sapi.SaveNexus(Filename=processed_filename, + InputWorkspace=self._samSqwWs) # additional output if self.getProperty("OutputSusceptibility").value: - temperature = mtd[self._samSqwWs].getRun().getProperty(TEMPERATURE_SENSOR).getStatistics().mean + temperature = mtd[self._samSqwWs].getRun().\ + getProperty(TEMPERATURE_SENSOR).getStatistics().mean samXqsWs = self._samSqwWs.replace("sqw", "Xqw") sapi.ApplyDetailedBalance(InputWorkspace=self._samSqwWs, - OutputWorkspace=samXqsWs, Temperature=str(temperature)) - sapi.ConvertUnits(InputWorkspace=samXqsWs, OutputWorkspace=samXqsWs, - Target="DeltaE_inFrequency", Emode="Indirect") + OutputWorkspace=samXqsWs, + Temperature=str(temperature)) + sapi.ConvertUnits(InputWorkspace=samXqsWs, + OutputWorkspace=samXqsWs, + Target="DeltaE_inFrequency", + Emode="Indirect") susceptibility_filename = processed_filename.replace("sqw", "Xqw") - sapi.SaveNexus(Filename=susceptibility_filename, InputWorkspace=samXqsWs) + sapi.SaveNexus(Filename=susceptibility_filename, + InputWorkspace=samXqsWs) if not self._debugMode: sapi.DeleteWorkspace("BASIS_MASK") # delete the mask if self._doNorm and bool(norm_runs): sapi.DeleteWorkspace("BASIS_NORM_MASK") # delete vanadium mask sapi.DeleteWorkspace(self._normWs) # Delete vanadium S(Q) + if self._normalizationType == "by Q slice": + sapi.DeleteWorkspace(normWs) # Delete vanadium events file def _getRuns(self, rlist, doIndiv=True): """ Create sets of run numbers for analysis. A semicolon indicates a separate group of runs to be processed together. @param rlist: string containing all the run numbers to be reduced. - @return if _doIndiv is False, return a list of IntArrayProperty objects. - Each item is a pseudolist containing a set of runs to be reduced together. - if _doIndiv is True, return a list of strings, each string is a run number. + @return if _doIndiv is False, return a list of IntArrayProperty + objects. Each item is a pseudolist containing a set of runs to + be reduced together. if _doIndiv is True, return a list of strings, + each string is a run number. """ run_list = [] - # ";" separates the runs into substrings. Each substring represents a set of runs + # ";" separates the runs into substrings. Each substring + # represents a set of runs rlvals = rlist.split(';') for rlval in rlvals: iap = IntArrayProperty("", rlval) # split the substring @@ -366,35 +411,58 @@ class BASISReduction(PythonAlgorithm): kwargs = {"BankName": "bank2"} # 311 analyzers only in bank2 else: kwargs = {} - sapi.LoadEventNexus(Filename=run_file, OutputWorkspace=ws_name, **kwargs) + sapi.LoadEventNexus(Filename=run_file, + OutputWorkspace=ws_name, **kwargs) if not self._noMonNorm: - sapi.LoadNexusMonitors(Filename=run_file, OutputWorkspace=mon_ws_name) + sapi.LoadNexusMonitors(Filename=run_file, + OutputWorkspace=mon_ws_name) if sam_ws != ws_name: - sapi.Plus(LHSWorkspace=sam_ws, RHSWorkspace=ws_name, OutputWorkspace=sam_ws) + sapi.Plus(LHSWorkspace=sam_ws, + RHSWorkspace=ws_name, + OutputWorkspace=sam_ws) sapi.DeleteWorkspace(ws_name) if mon_ws != mon_ws_name and not self._noMonNorm: - sapi.Plus(LHSWorkspace=mon_ws, RHSWorkspace=mon_ws_name, OutputWorkspace=mon_ws) + sapi.Plus(LHSWorkspace=mon_ws, + RHSWorkspace=mon_ws_name, + OutputWorkspace=mon_ws) sapi.DeleteWorkspace(mon_ws_name) def _calibData(self, sam_ws, mon_ws): - sapi.MaskDetectors(Workspace=sam_ws, DetectorList=self._dMask) # MaskedWorkspace='BASIS_MASK') - sapi.ModeratorTzeroLinear(InputWorkspace=sam_ws, OutputWorkspace=sam_ws) + sapi.MaskDetectors(Workspace=sam_ws, + DetectorList=self._dMask) + sapi.ModeratorTzeroLinear(InputWorkspace=sam_ws, + OutputWorkspace=sam_ws) sapi.LoadParameterFile(Workspace=sam_ws, - Filename=pjoin(DEFAULT_CONFIG_DIR, self._reflection["parameter_file"])) - sapi.ConvertUnits(InputWorkspace=sam_ws, OutputWorkspace=sam_ws, Target='Wavelength', EMode='Indirect') + Filename=pjoin(DEFAULT_CONFIG_DIR, + self._reflection["parameter_file"])) + sapi.ConvertUnits(InputWorkspace=sam_ws, + OutputWorkspace=sam_ws, + Target='Wavelength', + EMode='Indirect') if not self._noMonNorm: - sapi.ModeratorTzeroLinear(InputWorkspace=mon_ws, OutputWorkspace=mon_ws) - sapi.Rebin(InputWorkspace=mon_ws, OutputWorkspace=mon_ws, Params='10') - sapi.ConvertUnits(InputWorkspace=mon_ws, OutputWorkspace=mon_ws, Target='Wavelength') + sapi.ModeratorTzeroLinear(InputWorkspace=mon_ws, + OutputWorkspace=mon_ws) + sapi.Rebin(InputWorkspace=mon_ws, + OutputWorkspace=mon_ws, + Params='10') + sapi.ConvertUnits(InputWorkspace=mon_ws, + OutputWorkspace=mon_ws, + Target='Wavelength') sapi.OneMinusExponentialCor(InputWorkspace=mon_ws, OutputWorkspace=mon_ws, C='0.20749999999999999', C1='0.001276') - sapi.Scale(InputWorkspace=mon_ws, OutputWorkspace=mon_ws, Factor='1e-06') - sapi.RebinToWorkspace(WorkspaceToRebin=sam_ws, WorkspaceToMatch=mon_ws, OutputWorkspace=sam_ws) - sapi.Divide(LHSWorkspace=sam_ws, RHSWorkspace=mon_ws, OutputWorkspace=sam_ws) + sapi.Scale(InputWorkspace=mon_ws, + OutputWorkspace=mon_ws, + Factor='1e-06') + sapi.RebinToWorkspace(WorkspaceToRebin=sam_ws, + WorkspaceToMatch=mon_ws, + OutputWorkspace=sam_ws) + sapi.Divide(LHSWorkspace=sam_ws, + RHSWorkspace=mon_ws, + OutputWorkspace=sam_ws) def _sum_and_calibrate(self, run_set, extra_extension=""): """ @@ -420,9 +488,16 @@ class BASISReduction(PythonAlgorithm): @param isSample: discriminates between sample and vanadium @return: string name of S(Q,E) """ - sapi.ConvertUnits(InputWorkspace=wsName, OutputWorkspace=wsName, Target='DeltaE', EMode='Indirect') - sapi.CorrectKiKf(InputWorkspace=wsName, OutputWorkspace=wsName, EMode='Indirect') - sapi.Rebin(InputWorkspace=wsName, OutputWorkspace=wsName, Params=etRebins) + sapi.ConvertUnits(InputWorkspace=wsName, + OutputWorkspace=wsName, + Target='DeltaE', + EMode='Indirect') + sapi.CorrectKiKf(InputWorkspace=wsName, + OutputWorkspace=wsName, + EMode='Indirect') + sapi.Rebin(InputWorkspace=wsName, + OutputWorkspace=wsName, + Params=etRebins) if self._groupDetOpt != "None": if self._groupDetOpt == "Low-Resolution": grp_file = "BASIS_Grouping_LR.xml" @@ -432,33 +507,44 @@ class BASISReduction(PythonAlgorithm): # location to search paths if self._overrideMask: config.appendDataSearchDir(DEFAULT_MASK_GROUP_DIR) - sapi.GroupDetectors(InputWorkspace=wsName, OutputWorkspace=wsName, MapFile=grp_file, Behaviour="Sum") + sapi.GroupDetectors(InputWorkspace=wsName, + OutputWorkspace=wsName, + MapFile=grp_file, + Behaviour="Sum") # Output NXSPE file (must be done before transforming the # vertical axis to point data) if isSample and self._nsxpe_do: - #extension = '_divided_sqw.nxspe' if self._doNorm else '_sqw.nxspe' - extension = '_sqw.nxspe' + extension = '.nxspe' run = mtd[wsName].getRun() if run.hasProperty(self._nxspe_psi_angle_log): - psi_angle_logproperty = run.getProperty(self._nxspe_psi_angle_log) + psi_angle_logproperty = \ + run.getProperty(self._nxspe_psi_angle_log) psi_angle = np.average(psi_angle_logproperty.value) psi_angle += self._nxspe_offset nxspe_filename = wsName + extension sapi.SaveNXSPE(InputWorkspace=wsName, Filename=nxspe_filename, Efixed=self._reflection['default_energy'], - Psi=psi_angle, KiOverKfScaling=1) + Psi=psi_angle, + KiOverKfScaling=1) else: - error_message = 'Runs have no log entry named {}'.format(self._nxspe_psi_angle_log) + error_message = 'Runs have no log entry named {}'\ + .format(self._nxspe_psi_angle_log) self.log().error(error_message) - wsSqwName = wsName + '_divided_sqw' if isSample and self._doNorm else wsName + '_sqw' - sapi.SofQW3(InputWorkspace=wsName, QAxisBinning=self._qBins, EMode='Indirect', - EFixed=self._reflection["default_energy"], OutputWorkspace=wsSqwName) + wsSqwName = wsName + '_divided_sqw' \ + if isSample and self._doNorm else wsName + '_sqw' + sapi.SofQW3(InputWorkspace=wsName, + QAxisBinning=self._qBins, + EMode='Indirect', + EFixed=self._reflection["default_energy"], + OutputWorkspace=wsSqwName) # Rebin the vanadium within the elastic line if not isSample: - sapi.Rebin(InputWorkspace=wsSqwName, OutputWorkspace=wsSqwName, Params=self._reflection["vanadium_bins"]) + sapi.Rebin(InputWorkspace=wsSqwName, + OutputWorkspace=wsSqwName, + Params=self._reflection["vanadium_bins"]) return wsSqwName def _ScaleY(self, wsName): @@ -469,7 +555,10 @@ class BASISReduction(PythonAlgorithm): """ workspace = sapi.mtd[wsName] maximumYvalue = workspace.dataY(0).max() - sapi.Scale(InputWorkspace=wsName, OutputWorkspace=wsName, Factor=1./maximumYvalue, Operation="Multiply",) + sapi.Scale(InputWorkspace=wsName, + OutputWorkspace=wsName, + Factor=1./maximumYvalue, + Operation="Multiply") # Register algorithm with Mantid.