Newer
Older
Michael Whitty
committed
from mantidsimple import *
from mantidplot import *
Michael Whitty
committed
import IndirectEnergyConversion as IEC
Michael Whitty
committed
Michael Whitty
committed
import math
Michael Whitty
committed
import re
Michael Whitty
committed
def abscyl(inWS_n, outWS_n, efixed, sample, can):
Michael Whitty
committed
ConvertUnits(inWS_n, 'wavelength', 'Wavelength', 'Indirect', efixed)
CylinderAbsorption('wavelength', outWS_n, sample[0], sample[1], sample[2],
can[0],can[1], EMode='Indirect', EFixed=efixed, NumberOfSlices=can[2],
NumberOfAnnuli=can[3])
mantid.deleteWorkspace('wavelength')
Michael Whitty
committed
def absflat(inWS_n, outWS_n, efixed, sample, can):
Michael Whitty
committed
ConvertUnits(inWS_n, 'wlength', 'Wavelength', 'Indirect', efixed)
FlatPlateAbsorption('wlength', outWS_n, sample[0], sample[1], sample[2],
can[0], can[1], can[2], EMode='Indirect', EFixed=efixed,
ElementSize=can[3])
mantid.deleteWorkspace('wlength')
Michael Whitty
committed
Michael Whitty
committed
def absorption(input, mode, sample, can, Save=False, Verbose=False,
Michael Whitty
committed
Plot=False):
(direct, filename) = os.path.split(input)
(root, ext) = os.path.splitext(filename)
LoadNexusProcessed(input, root)
Michael Whitty
committed
ws = mtd[root]
det = ws.getDetector(0)
efixed = 0.0
try:
efixed = det.getNumberParameter('Efixed')[0]
except AttributeError: # detector group
ids = det.getDetectorIDs()
det = ws.getInstrument().getDetector(ids[0])
efixed = det.getNumberParameter('Efixed')[0]
Michael Whitty
committed
outWS_n = root[:-3] + 'abs'
if mode == 'Flat Plate':
absflat(root, outWS_n, efixed, sample, can)
if mode == 'Cylinder':
abscyl(root, outWS_n, efixed, sample, can)
if Save:
SaveNexus(outWS_n, outWS_n+'.nxs')
mantid.deleteWorkspace(root)
if Plot:
graph = plotSpectrum(outWS_n,0)
Michael Whitty
committed
def concatWSs(workspaces, unit, name):
dataX = []
dataY = []
dataE = []
for ws in workspaces:
readX = mtd[ws].readX(0)
readY = mtd[ws].readY(0)
readE = mtd[ws].readE(0)
for i in range(0, len(readX)):
dataX.append(readX[i])
for i in range(0, len(readY)):
dataY.append(readY[i])
for i in range(0, len(readE)):
dataE.append(readE[i])
CreateWorkspace(name, dataX, dataY, dataE, NSpec = len(workspaces),
UnitX=unit)
Michael Whitty
committed
def confitParsToWS(Table, Data, BackG='FixF'):
Michael Whitty
committed
dataX = []
ConvertSpectrumAxis(Data, 'inq', 'MomentumTransfer', 'Indirect')
Transpose('inq', 'inq')
readX = mtd['inq'].readX(0)
nBins = len(readX)
for i in range(0,nBins):
dataX.append(readX[i])
mtd.deleteWorkspace('inq')
xAxisVals = []
dataY = []
dataE = []
names = ''
ws = mtd[Table]
cName = ws.getColumnNames()
nSpec = ( ws.getColumnCount() - 1 ) / 2
for spec in range(0,nSpec):
yAxis = cName[(spec*2)+1]
if re.search('HWHM$', yAxis) or re.search('Height$', yAxis):
xAxisVals += dataX
if (len(names) > 0):
names += ","
names += yAxis
eAxis = cName[(spec*2)+2]
for row in range(0, ws.getRowCount()):
dataY.append(ws.getDouble(yAxis,row))
dataE.append(ws.getDouble(eAxis,row))
else:
nSpec -= 1
suffix = str(nSpec / 2) + 'L' + BackG
Michael Whitty
committed
outNm = Table + suffix
Michael Whitty
committed
CreateWorkspace(outNm, xAxisVals, dataY, dataE, nSpec,
Michael Whitty
committed
UnitX='MomentumTransfer', VerticalAxisUnit='Text',
VerticalAxisValues=names)
Michael Whitty
committed
return outNm
def confitPlotSeq(inputWS, plot):
nhist = mtd[inputWS].getNumberHistograms()
if ( plot == 'All' ):
plotSpectrum(inputWS, range(0, nhist))
return
plotSpecs = []
if ( plot == 'Intensity' ):
res = 'Height$'
elif ( plot == 'HWHM' ):
res = 'HWHM$'
for i in range(0,nhist):
title = mtd[inputWS].getAxis(1).label(i)
if re.search(res, title):
plotSpecs.append(i)
plotSpectrum(inputWS, plotSpecs)
Michael Whitty
committed
def confitSeq(inputWS, func, startX, endX, save, plot, bg):
Michael Whitty
committed
input = inputWS+',i0'
Michael Whitty
committed
for i in range(1, mtd[inputWS].getNumberHistograms()):
Michael Whitty
committed
input += ';'+inputWS+',i'+str(i)
Michael Whitty
committed
outNm = getWSprefix(inputWS) + '_conv_'
Michael Whitty
committed
PlotPeakByLogValue(input, outNm, func, StartX=startX, EndX=endX)
Michael Whitty
committed
wsname = confitParsToWS(outNm, inputWS, bg)
Michael Whitty
committed
if save:
SaveNexusProcessed(wsname, wsname+'.nxs')
Michael Whitty
committed
if plot != 'None':
confitPlotSeq(wsname, plot)
Michael Whitty
committed
def demon(rawfiles, first, last, instrument, Smooth=False, SumFiles=False,
grouping='Individual', cal='', CleanUp=True, Verbose=False,
Plot='None', Save=True, Real=False, Vanadium='', Monitor=True):
'''DEMON routine for diffraction reduction on indirect instruments (IRIS /
OSIRIS).'''
Michael Whitty
committed
# Get instrument workspace for gathering parameters and such
wsInst = mtd['__empty_' + instrument]
if wsInst is None:
IEC.loadInst(instrument)
wsInst = mtd['__empty_' + instrument]
# short name of instrument for saving etc
isn = ConfigService().facility().instrument(instrument).shortName().lower()
# parameters to do with monitor
if Monitor:
fmon, fdet = IEC.getFirstMonFirstDet('__empty_'+instrument)
try:
inst = wsInst.getInstrument()
area = inst.getNumberParameter('mon-area')[0]
thickness = inst.getNumberParameter('mon-thickness')[0]
except IndexError, message:
print message
sys.exit(message)
ws_mon_l = IEC.loadData(rawfiles, Sum=SumFiles, Suffix='_mon',
SpecMin=fmon+1, SpecMax=fmon+1)
Michael Whitty
committed
ws_det_l = IEC.loadData(rawfiles, Sum=SumFiles, SpecMin=first,
SpecMax=last)
Michael Whitty
committed
workspaces = []
for i in range(0, len(ws_det_l)):
det_ws = ws_det_l[i]
Michael Whitty
committed
# Get Run No
runNo = mtd[det_ws].getRun().getLogData("run_number").value
savefile = isn + runNo + '_diff'
if Monitor:
mon_ws = ws_mon_l[i]
# Get Monitor WS
IEC.timeRegime(monitor=mon_ws, detectors=det_ws, Smooth=Smooth)
IEC.monitorEfficiency(mon_ws, area, thickness)
IEC.normToMon(Factor=1e6, monitor=mon_ws, detectors=det_ws)
# Remove monitor workspace
DeleteWorkspace(mon_ws)
Michael Whitty
committed
if ( cal != '' ): # AlignDetectors and Group by .cal file
if Monitor:
ConvertUnits(det_ws, det_ws, 'TOF')
if ( mtd[det_ws].isDistribution() ):
ConvertFromDistribution(det_ws)
Michael Whitty
committed
AlignDetectors(det_ws, det_ws, cal)
DiffractionFocussing(det_ws, savefile, cal)
DeleteWorkspace(det_ws)
if ( Vanadium != '' ):
print "NotImplemented: divide by vanadium."
Michael Whitty
committed
else: ## Do it the old fashioned way
# Convert to dSpacing - need to AlignBins so we can group later
ConvertUnits(det_ws, det_ws, 'dSpacing', AlignBins=True)
IEC.groupData(grouping, savefile, detectors=det_ws)
Michael Whitty
committed
if Save:
SaveNexusProcessed(savefile, savefile+'.nxs')
Michael Whitty
committed
workspaces.append(savefile)
Michael Whitty
committed
if ( Plot != 'None' ):
Michael Whitty
committed
for demon in workspaces:
Michael Whitty
committed
if ( Plot == 'Contour' ):
importMatrixWorkspace(demon).plotGraph2D()
else:
nspec = mtd[demon].getNumberHistograms()
plotSpectrum(demon, range(0, nspec))
Michael Whitty
committed
return workspaces
Michael Whitty
committed
def elwin(inputFiles, eRange, Save=False, Verbose=False, Plot=False):
eq1 = [] # output workspaces with units in Q
eq2 = [] # output workspaces with units in Q^2
Michael Whitty
committed
for file in inputFiles:
(direct, filename) = os.path.split(file)
(root, ext) = os.path.splitext(filename)
LoadNexus(file, root)
Gigg, Martyn Anthony
committed
run = mtd[root].getRun().getLogData("run_number").value
Michael Whitty
committed
savefile = root[:3] + run + root[8:-3]
if ( len(eRange) == 4 ):
ElasticWindow(root, savefile+'eq1', savefile+'eq2',eRange[0],
eRange[1], eRange[2], eRange[3])
Michael Whitty
committed
elif ( len(eRange) == 2 ):
ElasticWindow(root, savefile+'eq1', savefile+'eq2',
eRange[0], eRange[1])
Michael Whitty
committed
if Save:
SaveNexusProcessed(savefile+'eq1', savefile+'eq1.nxs')
SaveNexusProcessed(savefile+'eq2', savefile+'eq2.nxs')
eq1.append(savefile+'eq1')
eq2.append(savefile+'eq2')
Michael Whitty
committed
mantid.deleteWorkspace(root)
if Plot:
Michael Whitty
committed
nBins = mtd[eq1[0]].getNumberBins()
lastXeq1 = mtd[eq1[0]].readX(0)[nBins-1]
graph1 = plotSpectrum(eq1, 0)
Michael Whitty
committed
layer = graph1.activeLayer()
layer.setScale(Layer.Bottom, 0.0, lastXeq1)
nBins = mtd[eq2[0]].getNumberBins()
lastXeq2 = mtd[eq2[0]].readX(0)[nBins-1]
graph2 = plotSpectrum(eq2, 0)
Michael Whitty
committed
layer = graph2.activeLayer()
layer.setScale(Layer.Bottom, 0.0, lastXeq2)
return eq1, eq2
Michael Whitty
committed
Michael Whitty
committed
def fury(sam_files, res_file, rebinParam, RES=True, Save=False, Verbose=False,
Plot=False):
outWSlist = []
# Process RES Data Only Once
LoadNexus(res_file, 'res_data') # RES
Rebin('res_data', 'res_data', rebinParam)
ExtractFFTSpectrum('res_data', 'res_fft', 2)
Integration('res_data', 'res_int')
Divide('res_fft', 'res_int', 'res')
for sam_file in sam_files:
(direct, filename) = os.path.split(sam_file)
(root, ext) = os.path.splitext(filename)
if (ext == '.nxs'):
LoadNexus(sam_file, 'sam_data') # SAMPLE
Rebin('sam_data', 'sam_data', rebinParam)
else: #input is workspace
Rebin(sam_file, 'sam_data', rebinParam)
Michael Whitty
committed
ExtractFFTSpectrum('sam_data', 'sam_fft', 2)
Integration('sam_data', 'sam_int')
Divide('sam_fft', 'sam_int', 'sam')
# Create save file name
savefile = getWSprefix('sam_data') + 'iqt'
Michael Whitty
committed
outWSlist.append(savefile)
Divide('sam', 'res', savefile)
#Cleanup Sample Files
mantid.deleteWorkspace('sam_data')
mantid.deleteWorkspace('sam_int')
mantid.deleteWorkspace('sam_fft')
mantid.deleteWorkspace('sam')
# Crop nonsense values off workspace
bin = int(math.ceil(mtd[savefile].getNumberBins()/ 2.0))
binV = mtd[savefile].dataX(0)[bin]
CropWorkspace(savefile, savefile, XMax=binV)
if Save:
SaveNexusProcessed(savefile, savefile + '.nxs')
# Clean Up RES files
mantid.deleteWorkspace('res_data')
mantid.deleteWorkspace('res_int')
mantid.deleteWorkspace('res_fft')
mantid.deleteWorkspace('res')
if Plot:
specrange = range(0,mtd[outWSlist[0]].getNumberHistograms())
plotFury(outWSlist, specrange)
return outWSlist
Michael Whitty
committed
Michael Whitty
committed
def furyfitCreateXAxis(inputWS):
Michael Whitty
committed
result = []
ws = mtd[inputWS]
nHist = ws.getNumberHistograms()
inst = ws.getInstrument()
samplePos = inst.getSample().getPos()
beamPos = samplePos - inst.getSource().getPos()
for i in range(0,nHist):
detector = ws.getDetector(i)
Michael Whitty
committed
try:
efixed = detector.getNumberParameter("Efixed")[0]
except AttributeError: # Detector Group
Michael Whitty
committed
ids = detector.getDetectorIDs()
Michael Whitty
committed
det = inst.getDetector(ids[0])
efixed = det.getNumberParameter("Efixed")[0]
Michael Whitty
committed
theta = detector.getTwoTheta(samplePos, beamPos) / 2
lamda = math.sqrt(81.787/efixed)
q = 4 * math.pi * math.sin(theta) / lamda
result.append(q)
return result
def furyfitParsToWS(Table, Data):
Michael Whitty
committed
dataX = furyfitCreateXAxis(Data)
dataY = []
dataE = []
names = ""
Michael Whitty
committed
xAxisVals = []
ws = mtd[Table]
cCount = ws.getColumnCount()
rCount = ws.getRowCount()
cName = ws.getColumnNames()
nSpec = ( cCount - 1 ) / 2
xAxis = cName[0]
Michael Whitty
committed
stretched = 0
for spec in range(0,nSpec):
yAxis = cName[(spec*2)+1]
Michael Whitty
committed
if ( re.search('Intensity$', yAxis) or re.search('Tau$', yAxis)
or re.search('Beta$', yAxis) ):
xAxisVals += dataX
if (len(names) > 0):
names += ","
names += yAxis
eAxis = cName[(spec*2)+2]
for row in range(0, rCount):
dataY.append(ws.getDouble(yAxis,row))
dataE.append(ws.getDouble(eAxis,row))
Michael Whitty
committed
if ( re.search('Beta$', yAxis) ): # need to know how many of curves
stretched += 1 # are stretched exponentials
else:
nSpec -= 1
Michael Whitty
committed
suffix = ''
nE = ( nSpec / 2 ) - stretched
if ( nE > 0 ):
suffix += str(nE) + 'E'
if ( stretched > 0 ):
suffix += str(stretched) + 'S'
wsname = Table + suffix
Michael Whitty
committed
CreateWorkspace(wsname, xAxisVals, dataY, dataE, nSpec,
Michael Whitty
committed
UnitX='MomentumTransfer', VerticalAxisUnit='Text',
VerticalAxisValues=names)
Michael Whitty
committed
return wsname
def furyfitPlotSeq(inputWS, plot):
nHist = mtd[inputWS].getNumberHistograms()
if ( plot == 'All' ):
plotSpectrum(inputWS, range(0, nHist))
return
plotSpecs = []
Michael Whitty
committed
if ( plot == 'Intensity' ):
res = 'Intensity$'
Michael Whitty
committed
if ( plot == 'Tau' ):
res = 'Tau$'
elif ( plot == 'Beta' ):
res = 'Beta$'
for i in range(0, nHist):
title = mtd[inputWS].getAxis(1).label(i)
if ( re.search(res, title) ):
plotSpecs.append(i)
plotSpectrum(inputWS, plotSpecs)
def furyfitSeq(inputWS, func, startx, endx, save, plot):
input = inputWS+',i0'
nHist = mtd[inputWS].getNumberHistograms()
for i in range(1,nHist):
input += ';'+inputWS+',i'+str(i)
Michael Whitty
committed
outNm = getWSprefix(inputWS) + 'fury_'
Michael Whitty
committed
PlotPeakByLogValue(input, outNm, func, StartX=startx, EndX=endx)
wsname = furyfitParsToWS(outNm, inputWS)
if save:
SaveNexusProcessed(wsname, wsname+'.nxs')
if ( plot != 'None' ):
furyfitPlotSeq(wsname, plot)
Michael Whitty
committed
Michael Whitty
committed
def getWSprefix(workspace):
'''Returns a string of the form '<ins><run>_<analyser><refl>_' on which
all of our other naming conventions are built.'''
if workspace == '':
return ''
Michael Whitty
committed
ws = mtd[workspace]
ins = ws.getInstrument().getName()
ins = ConfigService().facility().instrument(ins).shortName().lower()
run = ws.getRun().getLogData('run_number').value
analyser = ws.getInstrument().getStringParameter('analyser')[0]
reflection = ws.getInstrument().getStringParameter('reflection')[0]
prefix = ins + run + '_' + analyser + reflection + '_'
return prefix
Michael Whitty
committed
def msdfit(inputs, startX, endX, Save=False, Verbose=False, Plot=False):
output = []
Michael Whitty
committed
for file in inputs:
(direct, filename) = os.path.split(file)
(root, ext) = os.path.splitext(filename)
LoadNexusProcessed(file, root)
outWS_n = root[:-3] + 'msd'
fit_alg = Linear(root, outWS_n, WorkspaceIndex=0, StartX=startX,
EndX=endX)
output.append(outWS_n)
Michael Whitty
committed
A0 = fit_alg.getPropertyValue("FitIntercept")
A1 = fit_alg.getPropertyValue("FitSlope")
title = 'Intercept: '+A0+' ; Slope: '+A1
if Plot:
graph=plotSpectrum([root,outWS_n],0, 1)
graph.activeLayer().setTitle(title)
if Save:
SaveNexusProcessed(outWS_n, outWS_n+'.nxs', Title=title)
return output
Michael Whitty
committed
def plotFury(inWS_n, spec):
Michael Whitty
committed
nbins = inWS.getNumberBins()
graph = plotSpectrum(inWS_n, spec)
layer = graph.activeLayer()
Michael Whitty
committed
layer.setScale(Layer.Left, 0, 1.0)
Michael Whitty
committed
Michael Whitty
committed
def plotInput(inputfiles,spectra=[]):
OneSpectra = False
Michael Whitty
committed
if len(spectra) != 2:
Michael Whitty
committed
spectra = [spectra[0], spectra[0]]
OneSpectra = True
Michael Whitty
committed
workspaces = []
for file in inputfiles:
(direct, filename) = os.path.split(file)
(root, ext) = os.path.splitext(filename)
Michael Whitty
committed
LoadNexusProcessed(file, root)
if not OneSpectra:
GroupDetectors(root, root,
DetectorList=range(spectra[0],spectra[1]+1) )
Michael Whitty
committed
workspaces.append(root)
if len(workspaces) > 0:
graph = plotSpectrum(workspaces,0)
Michael Whitty
committed
layer = graph.activeLayer().setTitle(", ".join(workspaces))