Commit 86c9a799 authored by Steve Williams's avatar Steve Williams
Browse files

Fixes known problems with the Excitations interface. This GUI front end now...

Fixes known problems with the Excitations interface. This GUI front end now calls generates a call to an existing Python function, which makes changing the behaviour and scripting much easier. Numerical difference with libISIS seem to be fixed. refs #1137 and #841
parent c3eab224
......@@ -62,6 +62,7 @@ void GetEi::init()
"An approximate value for the typical incident energy, energy of\n"
"neutrons leaving the source (meV)");
declareProperty("IncidentEnergy", -1.0, Direction::Output);
declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
m_fracCompl = 0.0;
}
......@@ -86,26 +87,28 @@ void GetEi::exec()
// we're assuming that the time units for the X-values in the workspace are micro-seconds
const double peakLoc0 = 1e6*timeToFly(dist2moni0, E_est);
// write a lot of stuff to the log at user level as it is very possible for fit routines not to the expected thing
g_log.information() << "Based on the user selected energy the first peak will be searched for at TOF " << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%" << std::endl;
g_log.information() << "Based on the user selected energy the first peak will be searched for at TOF " << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%\n";
const double peakLoc1 = 1e6*timeToFly(dist2moni1, E_est);
g_log.information() << "Based on the user selected energy the second peak will be searched for at TOF " << peakLoc1 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%" << std::endl;
g_log.information() << "Based on the user selected energy the second peak will be searched for at TOF " << peakLoc1 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%\n";
// get the histograms created by the monitors
std::vector<int> indexes = getMonitorSpecIndexs(inWS, mon1Spec, mon2Spec);
g_log.information() << "Looking for a peak in the first monitor spectrum, spectra index " << indexes[0] << std::endl;
double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
g_log.information() << "The first peak has been found at TOF = " << t_monitor0 << std::endl;
g_log.notice() << "The first peak has been found at TOF = " << t_monitor0 << " microseconds\n";
setProperty("FirstMonitorPeak", t_monitor0);
g_log.information() << "Looking for a peak in the second monitor spectrum, spectra index " << indexes[1] << std::endl;
double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << std::endl;
g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << " microseconds\n";
// assumes that the source and the both mintors lie on one straight line, the 1e-6 converts microseconds to seconds as the mean speed needs to be in m/s
double meanSpeed = (dist2moni1 - dist2moni0)/(1e-6*(t_monitor1 - t_monitor0));
// uses 0.5mv^2 to get the kinetic energy in joules which we then convert to meV
double E_i = neutron_E_At(meanSpeed)/PhysicalConstants::meV;
g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV" << " (your estimate was " << E_est << " meV)" << std::endl;
g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV" << " (your estimate was " << E_est << " meV)\n";
setProperty("IncidentEnergy", E_i);
}
......@@ -126,7 +129,7 @@ void GetEi::getGeometry(DataObjects::Workspace2D_const_sptr WS, int mon0Spec, in
std::vector<int> dets = WS->spectraMap().getDetectors(mon0Spec);
if ( dets.size() != 1 )
{
g_log.error() << "The detector for spectrum number " << mon0Spec << " was either not found or is a group, grouped monitors are not supported by this algorithm" << std::endl;
g_log.error() << "The detector for spectrum number " << mon0Spec << " was either not found or is a group, grouped monitors are not supported by this algorithm\n";
g_log.error() << "Error retrieving data for the first monitor" << std::endl;
throw std::bad_cast();
}
......@@ -137,8 +140,8 @@ void GetEi::getGeometry(DataObjects::Workspace2D_const_sptr WS, int mon0Spec, in
dets = WS->spectraMap().getDetectors(mon1Spec);
if ( dets.size() != 1 )
{
g_log.error() << "The detector for spectrum number " << mon1Spec << " was either not found or is a group, grouped monitors are not supported by this algorithm" << std::endl;
g_log.error() << "Error retrieving data for the second monitor" << std::endl;
g_log.error() << "The detector for spectrum number " << mon1Spec << " was either not found or is a group, grouped monitors are not supported by this algorithm\n";
g_log.error() << "Error retrieving data for the second monitor\n";
throw std::bad_cast();
}
det = WS->getInstrument()->getDetector(dets[0]);
......@@ -307,7 +310,7 @@ void GetEi::getPeakEstimates(double &height, int &centreInd, double &background)
throw std::invalid_argument("No peak was found or its height is less than the threshold of " + boost::lexical_cast<std::string>(PEAK_THRESH_H) + " times the mean background, was the energy estimate (" + getPropertyValue("EnergyEstimate") + " meV) close enough?");
}
g_log.debug() << "Peak position is the bin that has the maximum Y value in the monitor spectrum, which is at TOF " << (m_tempWS->readX(0)[centreInd]+m_tempWS->readX(0)[centreInd+1])/2 << " (peak height " << height << " counts/microsecond)" << std::endl;
g_log.debug() << "Peak position is the bin that has the maximum Y value in the monitor spectrum, which is at TOF " << (m_tempWS->readX(0)[centreInd]+m_tempWS->readX(0)[centreInd+1])/2 << " (peak height " << height << " counts/microsecond)\n";
}
/** Estimates the closest time, looking either or back, when the number of counts is
......@@ -339,7 +342,7 @@ double GetEi::findHalfLoc(MantidVec::size_type startInd, const double height, co
if ( std::abs(static_cast<int>(endInd - startInd)) < PEAK_THRESH_W )
{// we didn't find a significant peak
g_log.error() << "Likely precision problem or error, one half height distance is less than the threshold number of bins from the central peak: " << std::abs(static_cast<int>(endInd - startInd)) << "<" << PEAK_THRESH_W << ". Check the monitor peak" << std::endl;
g_log.error() << "Likely precision problem or error, one half height distance is less than the threshold number of bins from the central peak: " << std::abs(static_cast<int>(endInd - startInd)) << "<" << PEAK_THRESH_W << ". Check the monitor peak\n";
}
// we have a peak in range, do an area check to see if the peak has any significance
double hOverN = (height-noise)/noise;
......@@ -360,7 +363,7 @@ double GetEi::findHalfLoc(MantidVec::size_type startInd, const double height, co
halfTime -= deltaY/gradient;
}
g_log.debug() << "One half height point found at TOF = " << halfTime << " microseconds" << std::endl;
g_log.debug() << "One half height point found at TOF = " << halfTime << " microseconds\n";
return halfTime;
}
/** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
......
......@@ -72,7 +72,7 @@ void GroupDetectors2::exec()
// Bin boundaries need to be the same, so do the full check on whether they actually are
if (!API::WorkspaceHelpers::commonBoundaries(inputWS))
{
g_log.error() << "Can only group if the histograms have common bin boundaries" << std::endl;
g_log.error() << "Can only group if the histograms have common bin boundaries\n";
throw std::invalid_argument("Can only group if the histograms have common bin boundaries");
}
progress( m_FracCompl = CHECKBINS );
......@@ -100,8 +100,7 @@ void GroupDetectors2::exec()
DataObjects::Workspace2D_sptr outputWS = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create(inputWS,
m_GroupSpecInds.size()+ numUnGrouped, inputWS->readX(0).size(),
inputWS->blocksize())
);
inputWS->blocksize()) );
// prepare to move the requested histograms into groups, first estimate how long for progress reporting. +1 in the demonator gets rid of any divide by zero risk
double prog4Copy=( (1.0 - m_FracCompl)/(numInHists-unGroupedSet.size()+1) )*
......@@ -116,7 +115,7 @@ void GroupDetectors2::exec()
moveOthers(unGroupedSet, inputWS, outputWS, outIndex);
}
g_log.information() << name() << " algorithm has finished" << std::endl;
g_log.information() << name() << " algorithm has finished\n";
setProperty("OutputWorkspace",outputWS);
}
......@@ -157,7 +156,7 @@ void GroupDetectors2::getGroups(DataObjects::Workspace2D_const_sptr workspace,
{
WorkspaceHelpers::getIndicesFromSpectra(workspace,
spectraList, m_GroupSpecInds[0]);
g_log.debug() << "Converted " << spectraList.size() << " spectra numbers into spectra indexes to be combined" << std::endl;
g_log.debug() << "Converted " << spectraList.size() << " spectra numbers into spectra indices to be combined\n";
}
else if ( ! detectorList.empty() )
{// we are going to group on the basis of detector IDs, convert from detectors to spectra numbers
......@@ -165,17 +164,17 @@ void GroupDetectors2::getGroups(DataObjects::Workspace2D_const_sptr workspace,
//then from spectra numbers to indices
WorkspaceHelpers::getIndicesFromSpectra(
workspace, mySpectraList, m_GroupSpecInds[0]);
g_log.debug() << "Found " << m_GroupSpecInds[0].size() << " spectra indexes from the list of " << detectorList.size() << " detectors" << std::endl;
g_log.debug() << "Found " << m_GroupSpecInds[0].size() << " spectra indices from the list of " << detectorList.size() << " detectors\n";
}
else if ( ! indexList.empty() )
{
m_GroupSpecInds[0] = indexList;
g_log.debug() << "Read in " << m_GroupSpecInds[0].size() << " spectra indexs to be combined" << std::endl;
g_log.debug() << "Read in " << m_GroupSpecInds[0].size() << " spectra indices to be combined\n";
}
if ( m_GroupSpecInds[0].empty() )
{
g_log.information() << name() << ": File, WorkspaceIndexList, SpectraList, and DetectorList properties are all empty" << std::endl;
g_log.information() << name() << ": File, WorkspaceIndexList, SpectraList, and DetectorList properties are all empty\n";
throw std::invalid_argument("All list properties are empty, nothing to group");
}
......@@ -186,7 +185,7 @@ void GroupDetectors2::getGroups(DataObjects::Workspace2D_const_sptr workspace,
{
unUsedSpec[m_GroupSpecInds[0][i]] = USED;
}
else g_log.warning() << "Duplicate index found" << std::endl;
else g_log.warning() << "Duplicate index found\n";
}
}
/** Read the spectra numbers in from the input file (the file format is in the
......@@ -208,9 +207,9 @@ void GroupDetectors2::processFile(std::string fname,
// for error reporting keep a count of where we are reading in the file
int lineNum = 1;
if ( File.rdstate() & std::ios::failbit )
if (File.fail())
{
g_log.debug() << " file state failbit set after reading attempt" << std::endl;
g_log.debug() << " file state failbit set after read attempt\n";
throw Exception::FileError("Couldn't read file", fname);
}
g_log.debug() << " success opening input file " << fname << std::endl;
......@@ -239,7 +238,7 @@ void GroupDetectors2::processFile(std::string fname,
if ( m_GroupSpecInds.size() != static_cast<size_t>(totalNumberOfGroups) )
{
g_log.warning() << "The input file header states there are " << totalNumberOfGroups << " but the file contains " << m_GroupSpecInds.size() << " groups" << std::endl;
g_log.warning() << "The input file header states there are " << totalNumberOfGroups << " but the file contains " << m_GroupSpecInds.size() << " groups\n";
}
}
// add some more info to the error messages, including the line number, to help users correct their files. These problems should cause the algorithm to stop
......@@ -248,8 +247,10 @@ void GroupDetectors2::processFile(std::string fname,
g_log.debug() << "Exception thrown: " << e.what() << std::endl;
File.close();
std::string error(e.what() + std::string(" near line number ") + boost::lexical_cast<std::string>(lineNum));
if ( File.rdstate() & std::ios::failbit )
if (File.fail())
{
error = "Input output error while reading file ";
}
throw Exception::FileError(error, fname);
}
catch (boost::bad_lexical_cast &e)
......@@ -257,12 +258,14 @@ void GroupDetectors2::processFile(std::string fname,
g_log.debug() << "Exception thrown: " << e.what() << std::endl;
File.close();
std::string error(std::string("Problem reading integer value \"") + e.what() + std::string("\" near line number ") + boost::lexical_cast<std::string>(lineNum));
if ( File.rdstate() & std::ios::failbit )
if (File.fail())
{
error = "Input output error while reading file ";
}
throw Exception::FileError(error, fname);
}
File.close();
g_log.debug() << "Closed file " << fname << " after reading in " << m_GroupSpecInds.size() << " groups" << std::endl;
g_log.debug() << "Closed file " << fname << " after reading in " << m_GroupSpecInds.size() << " groups\n";
m_FracCompl += fileReadProg( m_GroupSpecInds.size(), specs2index.size() );
return;
}
......@@ -402,7 +405,7 @@ void GroupDetectors2::readSpectraIndexes(std::string line, std::map<int,int> &sp
}
else
{// the spectra was already included in a group
g_log.warning() << "Duplicate spectra number " << specNumbers[i] << " ignored in input file" << std::endl;
g_log.warning() << "Duplicate spectra number " << specNumbers[i] << " ignored in input file\n";
}
}
}
......@@ -437,7 +440,7 @@ int GroupDetectors2::formGroups( DataObjects::Workspace2D_const_sptr inputWS, Da
// Get a reference to the spectra map on the output workspace
API::SpectraDetectorMap &specDetecMap = outputWS->mutableSpectraMap();
g_log.debug() << name() << ": Preparing to group spectra into " << m_GroupSpecInds.size() << " groups" << std::endl;
g_log.debug() << name() << ": Preparing to group spectra into " << m_GroupSpecInds.size() << " groups\n";
// where we are copying spectra to, we start copying to the start of the output workspace
int outIndex = 0;
......@@ -486,7 +489,7 @@ int GroupDetectors2::formGroups( DataObjects::Workspace2D_const_sptr inputWS, Da
}
outIndex ++;
}
g_log.debug() << name() << " created " << outIndex << " new grouped spectra" << std::endl;
g_log.debug() << name() << " created " << outIndex << " new grouped spectra\n";
return outIndex;
}
/** Only to be used if the KeepUnGrouped property is true, moves the spectra that were not selected
......@@ -526,7 +529,7 @@ void GroupDetectors2::moveOthers(const std::set<int> &unGroupedSet, DataObjects:
interruption_point();
}
}
g_log.debug() << name() << " copied " << unGroupedSet.size()-1 << " ungrouped spectra" << std::endl;
g_log.debug() << name() << " copied " << unGroupedSet.size()-1 << " ungrouped spectra\n";
}
//RangeHelper
/** Expands any ranges in the input string, eg. "1 3-5 4" -> "1 3 4 5 4"
......
......@@ -2,13 +2,11 @@ from mantidsimple import *
# returns a string with is comma separated list of the elements in the tuple, array or comma separated string!
def listToString(list):
stringIt = str(list)
stringIt = str(list).strip() #remove any white space from the front and end
if stringIt[0] == '[' or stringIt[0] == '(' :
# we get here if we were passed an array or tuple
lastInd = len(stringIt)-1
lastInd = len(stringIt)-1 # we get here if we were passed an array or tuple
if stringIt[lastInd] == ']' or stringIt[lastInd] == ')':
# only need to knock the brackets off
stringIt = stringIt[1:(lastInd-1)]
stringIt = stringIt[1:(lastInd)] # only need to knock the brackets off
return stringIt
def stringToList(commaSeparated):
......@@ -20,24 +18,25 @@ def stringToList(commaSeparated):
return theList
#sum all the workspaces, when the workspaces are not summed single input files are specified in this file and the final Python script is made of many copies of this file
def sumWorkspaces(total, runNumbers):
if len(runNumbers) > 1:
for toAdd in runNumbers[ 1 : ] :
#workspaces are overwritten to save memory
instru.loadRunNumber(toAdd, conv.tempWS)
Plus(total, conv.tempWS, total)
mantid.deleteWorkspace(conv.tempWS)
def sumWorkspaces(total, instrument, runNumbers):
if len(runNumbers) > 1:
tempWS = 'CommonFuncs_sumWorkspaces_tempory'
for toAdd in runNumbers[ 1 : ] :
instrument.loadRun(toAdd, tempWS) #workspaces are overwriten to save memory
Plus(total, tempWS, total)
mantid.deleteWorkspace(tempWS) #from www.mantidproject.org/MantidPyFramework
#-- Functions to do with input files
# uses the extension to decide whether use LoadNexus or LoadRaw
def LoadNexRaw(filename, workspace):
# this removes everything after the last, partition always returns three strings
filename = filename.strip() #remove any white space from the front or end of the name
# this removes everything after the last '.' (rpartition always returns three strings)
extension = filename.rpartition('.')[2]
if (extension == 'nxs') | (extension == 'NXS') :
if (extension.lower() == 'nxs') :
#return the first property from the algorithm, which for LoadNexus is the output workspace
return LoadNexus(filename, workspace)[0]
if (extension == 'raw') | (extension == 'RAW') :
return LoadRaw(filename, workspace, LoadLogFiles=0)[0]
if (extension.lower() == 'raw') :
return LoadRaw(filename, workspace)[0]
#we shouldn't get to here, the function should have returned by now
raise Exception("Could not find a load function for file "+filename+", *.raw and *.nxs accepted")
......@@ -52,7 +51,7 @@ def getFileName(instrumentPref, runNumber):
#only raw files are supported at the moment
return instrumentPref + runNumber + '.raw'
def loadRunNumber(instru, runNum, workspace):
def loadRun(instru, runNum, workspace):
return LoadNexRaw(getFileName(instru, runNum), workspace)
def loadMask(MaskFilename):
......@@ -91,10 +90,10 @@ class defaults:
def getFileName(self, runNumber):
try :
number = int(str(runNumber), 10)
except TypeError :
# means we weren't passed a number assume it is a valid file name and return it unprocessed
return runNumber
except Exception :
# means we weren't passed a number assume it is a valid file name and return it unprocessed
return runNumber
return getFileName(self.instrument_pref, number)
def loadRunNumber(self, runNum, workspace):
return LoadNexRaw(self.getFileName(runNum), workspace)
\ No newline at end of file
def loadRun(self, runNum, workspace):
return LoadNexRaw(self.getFileName(runNum), workspace)
from mantidsimple import *
import CommonFunctions as common
# these are the defaults for different instruments
import defaults
import ExcitDefaults
import math
##Estimate the location of the monitor peak, used to help with normalise to monitor
def getPeakTime(Ei, monSpec, WS):
source = WS.getInstrument().getSource();
#??STEVES?? replace MARI specific code below
# dets = WS.spectraMap().getDetectors(monSpec);
# if len(dets) != 1 : raise Exception('Found a grouped monitor, this is not supported')
def NormaliseTo(reference, WS, fromTOF):
if (reference == 'monitor-monitor 1') :
# ??STEVES read this in from defaults
min = 1000-fromTOF
max = 2000-fromTOF
NormaliseToMonitor( InputWorkspace=WS, OutputWorkspace=WS, MonitorSpectrum=1, IntegrationRangeMin=min, IntegrationRangeMax=max)
elif reference == 'monitor-monitor 2' :
PEAKWIDTH = 4.0/100
# the origin of time was set to the location of the peak so to get it's area integrate around the origin, PEAKWIDTH fraction of total time of flight
min = (-PEAKWIDTH/2)*20000
max = (PEAKWIDTH/2)*20000
#!! MARI Specific code!
det = WS.getDetector(\
\
\
monSpec-1\
\
\
\
);
distance = det.getDistance(source)
# convert distance to time, knowingthe kinetic energy:
# E_KE = mv^2/2, s = vt
# t = s/v, v = sqrt(2*E_KE/m)
# t = s/sqrt(2*E_KE/m)
# convert E_KE to joules kg m^2 s^-2
NeutronMass = 1.674927211e-27 #kg
meV = 1.602176487e-22 #joules
Ei *= meV;
# return the calculated time in micro-seconds
return (1e6*distance)/math.sqrt(2*Ei/NeutronMass)
def NormaliseTo(reference, WS, energy=0):
if (reference == 'monitor-monitor 1') :
if (energy == 0) : raise Exception('A non-zero energy must be supplied for normalisation by monitor')
# set the region where the monitor is
min = getPeakTime(energy, 2, WS)*0.96
max = getPeakTime(energy, 2, WS)*1.04
NormaliseToMonitor( InputWorkspace=WS, OutputWorkspace=WS, MonitorSpectrum=2, IntegrationRangeMin=min, IntegrationRangeMax=max)
elif reference == 'monitor-monitor 2' :
if (energy == 0) : raise Exception('A non-zero energy must be supplied for normalisation by monitor')
# set the region where the monitor is
min = getPeakTime(energy, 3, WS)*0.96
max = getPeakTime(energy, 3, WS)*1.04
NormaliseToMonitor( InputWorkspace=WS, OutputWorkspace=WS, MonitorSpectrum=3, IntegrationRangeMin=min, IntegrationRangeMax=max)
elif reference == 'protons (uAh)' :
NormaliseByCurrent( InputWorkspace=WS, OutputWorkspace=WS )
elif reference != 'no normalization' :
raise Exception('Normalisation scheme ' + reference + ' not found. It must be one of monitor, current, peak or none')
def NormaliseToWhiteBeam(WBRun, toNorm, mapFile, detMask, rebinString) :
def NormaliseToWhiteBeam(WBRun, toNorm,mapFile, detMask, rebinString, prevNorm, fromTOF) :
theNorm = "_ETrans_norm_tempory_WS"
try:
common.LoadNexRaw(WBRun, theNorm)
LoadDetectorInfo(theNorm, WBRun)
NormaliseTo(prevNorm, mtd[theNorm], fromTOF)
ConvertUnits(theNorm, theNorm, "Energy", AlignBins=0)
......@@ -82,43 +49,15 @@ def NormaliseToWhiteBeam(WBRun, toNorm, mapFile, detMask, rebinString) :
finally:
mantid.deleteWorkspace(theNorm)
def NormaliseToWhiteBeamAndLoadMask(WBRun, toNorm, mapFile, rebinString, DetectorMask) :
theNorm = "_ETrans_norm_tempory_WS"
try:
common.LoadNexRaw(WBRun, theNorm)
LoadDetectorInfo(theNorm, WBRun)
ConvertUnits(theNorm, theNorm, "Energy", AlignBins = 0)
#this both integrates the workspace into one bin spectra and sets up common bin boundaries for all spectra
Rebin(theNorm, theNorm, rebinString)
# shouldn't we do the correction? It affects things when the angles are different
# DetectorEfficiencyCor(theNorm, theNorm, IncidentE)
if DetectorMask != '' :
FDOL = FindDetectorsOutsideLimits(InputWorkspace=DetectorMask,OutputWorkspace='_ETrans_loading_bad_detector_WS',HighThreshold=10,LowThreshold=-1,OutputFile='')
badNums = FDOL.getPropertyValue('BadSpectraNums')
MaskDetectors(Workspace=theNorm, SpectraList=badNums)
mantid.deleteWorkspace('_ETrans_loading_bad_detector_WS')
if mapFile!= "" :
GroupDetectors( theNorm, theNorm, mapFile, KeepUngroupedSpectra=0)
Divide(toNorm, theNorm, toNorm)
#bad detectors were identified by the last command as infinities or NaN (0/0 = Not a Number)
ReplaceSpecialValues(toNorm, toNorm, NaNValue=-1e30, InfinityValue=-1e30)
finally:
mantid.deleteWorkspace(theNorm)
# returns the background range from the string range, which could be two numbers separated by a coma or could be empty and then the default values in instrument are used. Another possiblity is that background correction isn't going to be done, then we don't have to do anything
# throws if the string range is in the wrong format
def getBackRange(range, instrument) :
if range == 'noback' : ['noback', 'noback']
if range == 'noback' : return ['noback', 'noback']
# load the default time of flight values that delimit the background region
TOFLow = instrument.backgroundRange[0]
TOFHigh = instrument.backgroundRange[1]
if range != '' :
range = common.listToString(range)
strings = range.split(',')
TOFLow = float(strings[0])
TOFHigh = float(strings[1])
......@@ -131,17 +70,40 @@ def getNorm(instrument, normMethod):
###--Read in user set parameters or defaults values---------------
def retrieveSettings(instrum, back, runs):
instru = defaults.loadDefaults(instrum)
instru = ExcitDefaults.loadDefaults(instrum)
run_nums = common.listToString(runs).split(',')
try:
[TOFLow, TOFHigh] = getBackRange(back, instru)
except:
raise Exception('The background range must be specified as a string\n Two numbers separated by a coma to specify the range\n noback to disable background correction\n an empty string or absent for no background correction')
run_nums = common.listToString(runs).split(',')
return (instru, TOFLow, TOFHigh, run_nums)
def findPeaksOffSetTimeAndEi(workS, E_guess, getEi=True, detInfoFile=''):
if detInfoFile != '' :
LoadDetectorInfo(workS, detInfoFile) #for details see www.mantidproject.org/LoadDetectorInfo
runGetEi = str(getEi).lower() != 'fixei' and str(getEi).lower() != 'false'
if not runGetEi :
mtd.sendLogMessage('Getting peak times (the GetEi() energy will be overriden by '+str(E_guess)+' meV)')
#??STEVES?? MARI specfic code GetEiData.FirstMonitor ?
GetEiData = GetEi(workS, 2, 3, E_guess)
if runGetEi : Ei = GetEiData.getPropertyValue('IncidentEnergy')
else : Ei = E_guess
offSetTime=float(GetEiData.getPropertyValue('FirstMonitorPeak')) # set the origin of time to be when the neutrons arrive at the first GetEi monitor (often called monitor 2)
ChangeBinOffset(workS, workS, -offSetTime)
#??STEVES?? MARI specfic code GetEiData.FirstMonitor ?
p0=workS.getDetector(1).getPos() # set the source to be where the neutrons were at the origin of time (that first GetEi monitor)
mtd.sendLogMessage('Adjusting the TOF values to start when the peak was registered by the monitor (2) at the new origin')
src = workS.getInstrument().getSource().getName() #this name could be 'Moderator', 'undulator', etc.
MoveInstrumentComponent(workS, src, X=p0.getX() , Y=p0.getY() ,Z=p0.getZ(), RelativePosition=False)
return [Ei, offSetTime]
# a workspace name that is very long and unlikely to have been created by the user before this script is run, it would be replaced
tempWS = '_ETrans_loading_tempory_WS'
......@@ -149,32 +111,24 @@ tempWS = '_ETrans_loading_tempory_WS'
# Applies unit conversion, detector convcy and grouping correction to a
# raw file
###########################################
def mono_sample(instrum, runs, Ei, d_rebin, wbrf, wb_low_e=0, wb_high_e=1e8, getEi='true', back='', norma='', det_map='', det_mask='', scaling=1):
nameInOut = 'mono_sample_temporyWS'
def mono_sample(instrum, runs, Ei, d_rebin, wbrf, wb_low_e=0, wb_high_e=1e8, getEi=True, back='', norma='', det_map='', det_mask='', nameInOut = 'mono_sample_temporyWS') :
(instru, TOFLow, TOFHigh, run_nums) = retrieveSettings(instrum, back, runs)
#----Calculations start------------------
try:
#this new workspace is accessed by its name nameInOut = 'mono_sample_temporyWS' or pInOut (a pointer)
pInOut = instru.loadRunNumber(run_nums[0], nameInOut)
pInOut = instru.loadRun(run_nums[0], nameInOut) #this new workspace is accessed by its name nameInOut = 'mono_sample_temporyWS' or pInOut (a pointer)
if pInOut.isGroup() : raise Exception("Workspace groups are not supported here")
if getEi.lower() != 'fixei' and getEi.lower() != 'false' :
#?? 2 and 3 below are good for MARI, this needs to be changed in other instruments
GetEiData = GetEi(pInOut, 2, 3, Ei)
Ei = GetEiData.getPropertyValue('IncidentEnergy')
LoadDetectorInfo(pInOut, instru.getFileName(run_nums[0]))
common.sumWorkspaces(pInOut, run_nums)
if back != 'noback':
ConvertToDistribution(pInOut)
FlatBackground(pInOut, pInOut, TOFLow, TOFHigh, '', 'Mean')
ConvertFromDistribution(pInOut)
common.sumWorkspaces(pInOut, instru, run_nums) #the final workspace will take the X-values and instrument from the first workspace and so we don't need to rerun ChangeBinOffset(), MoveInstrumentComponent(), etc.
[Ei, offSet] = findPeaksOffSetTimeAndEi(pInOut, Ei, getEi, instru.getFileName(run_nums[0]))
if back != 'noback': #remove the count rate seen in the regions of the histograms defined as the background regions, if the user defined a region
ConvertToDistribution(pInOut) #deal correctly with changing bin widths
FlatBackground(pInOut, pInOut, TOFLow-offSet, TOFHigh-offSet, '', 'Mean') #the TOFs in the workspace were offset in a command above so we must take care referring to times now
ConvertFromDistribution(pInOut) #some algorithms later can't deal with distributions
# deals with normalize to monitor (e.g. norm = 'monitor-monitor 1'), current (if norm = 'protons (uAh)'), etc.
NormaliseTo(getNorm(instru, norma), pInOut, float(Ei))
NormaliseTo(getNorm(instru, norma), pInOut, offSet)
ConvertUnits(pInOut, pInOut, 'DeltaE', 'Direct', Ei, AlignBins=0)
......@@ -191,18 +145,15 @@ def mono_sample(instrum, runs, Ei, d_rebin, wbrf, wb_low_e=0, wb_high_e=1e8, get
ConvertToDistribution(pInOut)
# multiply by the user defined arbitary scaling factor, used because some plotting applications prefer numbers close to 1
if scaling != 1:
CreateSingleValuedWorkspace(tempWS, scaling)
Multiply(pInOut, tempWS, pInOut)
mantid.deleteWorkspace(tempWS)
CreateSingleValuedWorkspace(tempWS, 1.8182e8)
Multiply(pInOut, tempWS, pInOut)
mantid.deleteWorkspace(tempWS)
#replaces inifinities and error values with large numbers. Infinity values can be normally be avoided passing good energy values to ConvertUnits
ReplaceSpecialValues(pInOut, pInOut, 1e40, 1e40, 1e40, 1e40)
wb_rebin = str(wb_low_e) + ', ' + str(2*float(wb_high_e)) + ', ' + str(wb_high_e)
NormaliseToWhiteBeam(instru.getFileName(wbrf), pInOut, det_map, det_mask, wb_rebin)
#bad detectors were identified by the last command as infinities or NaN (0/0 = Not a Number)
ReplaceSpecialValues(pInOut, pInOut, NaNValue=-1e30, InfinityValue=-1e30)
NormaliseToWhiteBeam(instru.getFileName(wbrf), pInOut, det_map, det_mask, wb_rebin, norma, offSet)
return pInOut
......@@ -214,3 +165,7 @@ def mono_sample(instrum, runs, Ei, d_rebin, wbrf, wb_low_e=0, wb_high_e=1e8, get
print reason
#best to raise the exception if you know that all your Python environments can take that
#raise
#below is a quick test/example command to run this library, it must be commented or the library won't work!
#badSpectra = diagnose(instrum='MAR',wbrf='11060',wbrf2='11060',runs='15537',tiny=1e-10,huge=1e10,median_lo=0.1,median_hi=3.0,sv_sig=3.3,bmedians=5.0,zero='False', out_asc='',maskFile='')
#pRes = mono_sample('MAR', '15537', 82.541, '-10,0.1, 70', 11060, wb_low_e=20, wb_high_e=40, getEi='false', back='', norma='monitor-monitor 1', det_map='mari_res.map', det_mask=[])
\ No newline at end of file
......@@ -7,7 +7,7 @@ from os import remove
from mantidsimple import *
import CommonFunctions as common
# these are the defaults for different instruments
import defaults
import ExcitDefaults
OUT_WS_PREFIX = 'mask_'
......@@ -109,7 +109,7 @@ def single_white_van(wbrf, tiny, huge, median_hi, median_lo, erro