IndirectDataAnalysis.py

from mantidsimple import *
from mantidplot import *
from IndirectEnergyConversion import *
import math

def abscyl(inWS_n, outWS_n, efixed, sample, can):
    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')

def absflat(inWS_n, outWS_n, efixed, sample, can):
    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')

def absorption(input, mode, sample, can, efixed, Save=False, Verbose=False,
        Plot=False):
    (direct, filename) = os.path.split(input)
    (root, ext) = os.path.splitext(filename)
    LoadNexusProcessed(input, root)
    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)

def demon(rawFiles, first, last, Smooth=False, SumFiles=False, CleanUp=True,
        Verbose=False, Plot=False, Save=True):
    ws_list, ws_names = loadData(rawFiles, Sum=SumFiles)
    runNos = []
    workspaces = []
    (direct, filename) = os.path.split(rawFiles[0])
    (root, ext) = os.path.splitext(filename)
    try:
        inst = mtd[ws_names[0]].getInstrument()
        area = inst.getNumberParameter('mon-area')[0]
        thickness = inst.getNumberParameter('mon-thickness')[0]
    except IndexError:
        sys.exit('Monitor area and thickness (unt and zz) are not defined \
                in the Instrument Parameter File.')
    for i in range(0, len(ws_names)):
        # Get Monitor WS
        MonitorWS = timeRegime(inWS=ws_names[i], Smooth=Smooth)
        monitorEfficiency(MonitorWS, area, thickness)
        # Get Run no, crop file
        runNo = ws_list[i].getRun().getLogData("run_number").value()
        runNos.append(runNo)
        savefile = root[:3] + runNo + '_dem'
        CropWorkspace(ws_names[i], ws_names[i], StartWorkspaceIndex=(first-1),
                EndWorkspaceIndex = (last-1) )
        # Normalise to Monitor
        normalised = normToMon(inWS_n=ws_names[i], outWS_n=ws_names[i],
            monWS_n=MonitorWS)
        # Convert to dSpacing
        ConvertUnits(ws_names[i], savefile, 'dSpacing')
        workspaces.append(savefile)
        if CleanUp:
            mantid.deleteWorkspace(ws_names[i])
        if Save:
            SaveNexusProcessed(savefile, savefile+'.nxs')
    if Plot:
        for demon in workspaces:
            nspec = mantid.getMatrixWorkspace(demon).getNumberHistograms()
            plotSpectrum(demon, range(0, nspec))
    return workspaces, runNos

def elwin(inputFiles, eRange, efixed, Save=False, Verbose=False, Plot=False):
    outWS_list = []
    for file in inputFiles:
        (direct, filename) = os.path.split(file)
        (root, ext) = os.path.splitext(filename)
        LoadNexus(file, root)
        run = mtd[root].getRun().getLogData("run_number").value()
        savefile = root[:3] + run + '_elw'
        Elwin(root, savefile, eRange, efixed)
        if Save:
            SaveNexusProcessed(savefile, savefile+'.nxs')
        outWS_list.append(savefile)
        mantid.deleteWorkspace(root)
    if Plot:
        graph = plotSpectrum(outWS_list, 0)
    return outWS_list

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)
        LoadNexus(sam_file, 'sam_data') # SAMPLE
        Rebin('sam_data', 'sam_data', rebinParam)
        ExtractFFTSpectrum('sam_data', 'sam_fft', 2)
        Integration('sam_data', 'sam_int')
        Divide('sam_fft', 'sam_int', 'sam')
        # Create save file name.
        runNo = mtd['sam_data'].getRun().getLogData("run_number").value()
        savefile = root[:3] + runNo + '_iqt'
        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

def furyfitSeq(inputWS,func,startx,endx):
    input = inputWS+',i0'
    nHist = mtd[inputWS].getNumberHistograms()
    for i in range(1,nHist):
        input += ';'+inputWS+',i'+str(i)
    outNm = inputWS + '_fitParameters'
    PlotPeakByLogValue(input, outNm, func, StartX=startx, EndX=endx)
    procSeqParToWS(outNm)

def procSeqParToWS(inputWS):
    ws = mtd[inputWS]
    dataX = []
    dataY = []
    dataE = []
    cCount = ws.getColumnCount()
    rCount = ws.getRowCount()
    cName =  ws.getColumnNames()
    nSpec = ( cCount - 1 ) / 2
    xAxis = cName[0]
    for spec in range(0,nSpec):
        yAxis = cName[(spec*2)+1]
        eAxis = cName[(spec*2)+2]
        for row in range(0, rCount):
            dataX.append(ws.getDouble(xAxis,row))
            dataY.append(ws.getDouble(yAxis,row))
            dataE.append(ws.getDouble(eAxis,row))
    CreateWorkspace(inputWS+'_matrix', dataX, dataY, dataE, nSpec)

def msdfit(inputs, startX, endX, Save=False, Verbose=False, Plot=False):
    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)
        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)

def mut(inWS_n, deltaW, filename, efixed):
    file_handle = open(filename, 'w') # Open File
    sigs = 5.0 # sigs ?
    siga = 1.0 # siga ?
    we = math.sqrt(81.787/efixed) # 81.787 ?
    tempWS = '_tmp_indirect_mut_file_'
    ConvertUnits(inWS_n, tempWS, 'Wavelength', 'Indirect', efixed)
    xValues = mtd[tempWS].readX(0)
    nbins = len(xValues)
    xMin = xValues[0]
    xMax = xValues[nbins-1]
    nw = int(xMax/deltaW) - int(xMin/deltaW) + 1
    file_handle.write(str(nw)+' '+str(we)+' '+str(siga)+' '+str(sigs)+'\n')
    for i in range(0, nw):
        wavelength = (int(xMin/deltaW) + i) * deltaW
        sigt = sigs + siga*wavelength/1.8
        file_handle.write(str(wavelength) + ' ' + str(sigt) + '\n')
    file_handle.close()
    mantid.deleteWorkspace(tempWS)

def plotFury(inWS_n, spec):
    inWS = mantid.getMatrixWorkspace(inWS_n[0])
    nbins = inWS.getNumberBins()
    lastValueZero = False
    for i in range(0, nbins):
        if (inWS.readY(spec[0])[i] != inWS.readY(spec[0])[i]):
            if lastValueZero:
                xBoundary = inWS.readX(spec[0])[i]
                break
            else:
                lastValueZero = True
        else:
            lastValueZero = False
    graph = plotSpectrum(inWS_n, spec)
    layer = graph.activeLayer()
    layer.setScale(0, 0, 1.0)
    layer.setScale(2, 0, xBoundary)

def slice(inputfiles, calib, xrange, spec, Save=False, Verbose=False,
        Plot=False):
    outWSlist = []
    if  not ( ( len(xrange) == 2 ) or ( len(xrange) == 4 ) ):
        mantid.sendLogMessage('>> TOF Range must contain either 2 or 4 \
                numbers.')
        sys.exit(1)
    for file in inputfiles:
        (direct, filename) = os.path.split(file)
        (root, ext) = os.path.splitext(filename)
        if spec == [0, 0]:
            LoadRaw(file, root)
        else:
            LoadRaw(file, root, SpectrumMin=spec[0], SpectrumMax=spec[1])
        nhist = mtd[root].getNumberHistograms()
        if calib != '':
            useCalib(calib, inWS_n=root, outWS_n=root)
        run = mtd[root].getRun().getLogData("run_number").value()
        sfile = root[:3] + run + '_slt'
        if (len(xrange) == 2):
            Integration(root, sfile, xrange[0], xrange[1], 0, nhist-1)
        else:
            FlatBackground(root, sfile, StartX=xrange[2], EndX=xrange[3], 
                    Mode='Mean')
            Integration(sfile, sfile, xrange[0], xrange[1], 0, nhist-1)
        if Save:
            SaveNexusProcessed(sfile, sfile+'.nxs')
        outWSlist.append(sfile)
        mantid.deleteWorkspace(root)
    if Plot:
        graph = plotBin(outWSlist, 0)

def plotRaw(inputfiles,spectra=[]):
    if len(spectra) != 2:
        sys.exit(1)
    workspaces = []
    for file in inputfiles:
        (direct, filename) = os.path.split(file)
        (root, ext) = os.path.splitext(filename)
        LoadRaw(file, root, SpectrumMin=spectra[0], SpectrumMax = spectra[1])
        GroupDetectors(root,root,DetectorList=range(spectra[0],spectra[1]+1))
        workspaces.append(root)
    if len(workspaces) > 0:
        graph = plotSpectrum(workspaces,0)
        layer = graph.activeLayer().setTitle(", ".join(workspaces))

def plotInput(inputfiles,spectra=[]):
    if len(spectra) != 2:
        sys.exit(1)
    workspaces = []
    for file in inputfiles:
        (direct, filename) = os.path.split(file)
        (root, ext) = os.path.splitext(filename)
        LoadNexusProcessed(file, root)
        GroupDetectors(root,root,DetectorList=range(spectra[0],spectra[1]+1))
        workspaces.append(root)
    if len(workspaces) > 0:
        graph = plotSpectrum(workspaces,0)
        layer = graph.activeLayer().setTitle(", ".join(workspaces))