Ticket #8390: IndirectDataAnalysis.py

File IndirectDataAnalysis.py, 51.5 KB (added by Samuel Jackson, 7 years ago)

Updated version of the script

Line 
1from mantid.simpleapi import *
2from IndirectImport import import_mantidplot
3mp = import_mantidplot()
4from IndirectCommon import *
5from mantid import config, logger
6import math, re, os.path, numpy as np
7
8##############################################################################
9# Misc. Helper Functions
10##############################################################################
11
12def split(l, n):
13    #Yield successive n-sized chunks from l.
14    for i in xrange(0, len(l), n):
15        yield l[i:i+n]
16
17def segment(l, fromIndex, toIndex):
18    for i in xrange(fromIndex, toIndex + 1):
19        yield l[i]
20
21def trimData(nSpec, vals, min, max):
22    result = []
23    chunkSize = len(vals) / nSpec
24    assert min >= 0, 'trimData: min is less then zero'
25    assert max <= chunkSize - 1, 'trimData: max is greater than the number of spectra'
26    assert min <= max, 'trimData: min is greater than max'
27    chunks = split(vals,chunkSize)
28    for chunk in chunks:
29        seg = segment(chunk,min,max)
30        for val in seg:
31            result.append(val)
32    return result
33
34##############################################################################
35# ConvFit
36##############################################################################
37
38def getConvFitOption(ftype, bgd, Verbose):
39    if ftype[:5] == 'Delta':
40        delta = True
41        lor = ftype[5:6]
42    else:
43        delta = False
44        lor = ftype[:1]
45    options = [bgd, delta, int(lor)]
46    if Verbose:
47        logger.notice('Fit type : Delta = ' + str(options[1]) + ' ; Lorentzians = ' + str(options[2]))
48        logger.notice('Background type : ' + options[0])
49    return options
50
51##############################################################################
52
53def createConvFitFun(options, par, file):
54    bgd_fun = 'name=LinearBackground,A0='
55    if options[0] == 'FixF':
56        bgd_fun = bgd_fun +str(par[0])+',A1=0,ties=(A0='+str(par[0])+',A1=0.0)'
57    if options[0] == 'FitF':
58        bgd_fun = bgd_fun +str(par[0])+',A1=0,ties=(A1=0.0)'
59    if options[0] == 'FitL':
60        bgd_fun = bgd_fun +str(par[0])+',A1='+str(par[1])
61    if options[1]:
62        ip = 3
63    else:
64        ip = 2
65    pk_1 = '(composite=Convolution;name=Resolution, FileName="'+file+'"'
66    if  options[2] >= 1:
67        lor_fun = 'name=Lorentzian,Amplitude='+str(par[ip])+',PeakCentre='+str(par[ip+1])+',HWHM='+str(par[ip+2])
68    if options[2] == 2:
69        funcIndex = 1 if options[1] else 0
70        lor_2 = 'name=Lorentzian,Amplitude='+str(par[ip+3])+',PeakCentre='+str(par[ip+4])+',HWHM='+str(par[ip+5])
71        lor_fun = lor_fun +';'+ lor_2 +';ties=(f'+str(funcIndex)+'.PeakCentre=f'+str(funcIndex+1)+'.PeakCentre)'
72    if options[1]:
73        delta_fun = 'name=DeltaFunction,Height='+str(par[2])
74        lor_fun = delta_fun +';' + lor_fun
75    func = bgd_fun +';'+ pk_1 +';('+ lor_fun +'))'
76    return func
77
78##############################################################################
79
80def getConvFitResult(inputWS, resFile, outNm, ftype, bgd, specMin, specMax, Verbose):
81    options = getConvFitOption(ftype, bgd[:-2], Verbose)   
82    params = mtd[outNm+'_Parameters']
83    A0 = params.column(1)     #bgd A0 value
84    A1 = params.column(3)     #bgd A1 value
85    if options[1]:
86        ip = 7
87        D1 = params.column(5)      #delta value
88    else:
89        ip = 5
90    if options[2] >= 1:
91        H1 = params.column(ip)        #height1 value
92        C1 = params.column(ip+2)      #centre1 value
93        W1 = params.column(ip+4)      #width1 value
94    if options[2] == 2:
95        H2 = params.column(ip+6)        #height2 value
96        C2 = params.column(ip+8)      #centre2 value
97        W2 = params.column(ip+10)      #width2 value
98
99    for i in range(0,specMax-specMin):
100        paras = [A0[i], A1[i]]
101        if options[1]:
102            paras.append(D1[i])
103        if options[2] >= 1:
104            paras.append(H1[i])
105            paras.append(C1[i])
106            paras.append(W1[i])
107        if options[2] == 2:
108            paras.append(H2[i])
109            paras.append(C2[i])
110            paras.append(W2[i])
111        func = createConvFitFun(options, paras, resFile)
112        if Verbose:
113            logger.notice('Fit func : '+func)     
114        fitWS = outNm + '_Result_'
115        fout = fitWS + str(i)
116        Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i+specMin,Output=fout,MaxIterations=0)
117        unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
118        unitx.setLabel('Time' , 'ns')
119        RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
120        AddSampleLog(Workspace=fout, LogName="Fit Program", LogType="String", LogText='ConvFit')
121        AddSampleLog(Workspace=fout, LogName='Background', LogType='String', LogText=str(options[0]))
122        AddSampleLog(Workspace=fout, LogName='Delta', LogType='String', LogText=str(options[1]))
123        AddSampleLog(Workspace=fout, LogName='Lorentzians', LogType='String', LogText=str(options[2]))
124        DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
125        DeleteWorkspace(fitWS+str(i)+'_Parameters')
126        if i == 0:
127            group = fout
128        else:
129            group += ',' + fout
130    GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
131
132##############################################################################
133
134def confitParsToWS(Table, Data, specMin=0, specMax=-1):
135    if ( specMax == -1 ):
136        specMax = mtd[Data].getNumberHistograms() - 1
137    dataX = createQaxis(Data)
138    xAxisVals = []
139    xAxisTrimmed = []
140    dataY = []
141    dataE = []
142    names = ''
143    ws = mtd[Table]
144    cName =  ws.getColumnNames()
145    nSpec = ( ws.columnCount() - 1 ) / 2
146    for spec in range(0,nSpec):
147        yCol = (spec*2)+1
148        yAxis = cName[(spec*2)+1]
149        if re.search('HWHM$', yAxis) or re.search('Amplitude$', yAxis):
150            xAxisVals += dataX
151            if (len(names) > 0):
152                names += ","
153            names += yAxis
154            eCol = (spec*2)+2
155            eAxis = cName[(spec*2)+2]
156            for row in range(0, ws.rowCount()):
157                dataY.append(ws.cell(row,yCol))
158                dataE.append(ws.cell(row,eCol))
159        else:
160            nSpec -= 1
161    outNm = Table + "_Workspace"
162    xAxisTrimmed = trimData(nSpec, xAxisVals, specMin, specMax)
163    CreateWorkspace(OutputWorkspace=outNm, DataX=xAxisTrimmed, DataY=dataY, DataE=dataE, 
164        Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
165        VerticalAxisValues=names)
166    return outNm
167
168##############################################################################
169
170def confitPlotSeq(inputWS, Plot):
171    nhist = mtd[inputWS].getNumberHistograms()
172    if ( Plot == 'All' ):
173        mp.plotSpectrum(inputWS, range(0, nhist), True)
174        return   
175    plotSpecs = []
176    if ( Plot == 'Intensity' ):
177        res = 'Amplitude$'
178    elif ( Plot == 'HWHM' ):
179        res = 'HWHM$'
180    for i in range(0,nhist):
181        title = mtd[inputWS].getAxis(1).label(i)
182        if re.search(res, title):
183            plotSpecs.append(i)
184    mp.plotSpectrum(inputWS, plotSpecs, True)
185
186##############################################################################
187
188def confitSeq(inputWS, func, startX, endX, Save, Plot, ftype, bgd, specMin, specMax, Verbose):
189    StartTime('ConvFit')
190    workdir = config['defaultsave.directory']
191    elements = func.split('"')
192    resFile = elements[1] 
193    if Verbose:
194        logger.notice('Input files : '+str(inputWS)) 
195    input = inputWS+',i' + str(specMin)
196    if (specMax == -1):
197        specMax = mtd[inputWS].getNumberHistograms() - 1
198    for i in range(specMin + 1, specMax + 1):
199        input += ';'+inputWS+',i'+str(i)
200    (instr, run) = getInstrRun(inputWS)
201    run_name = instr + run
202    outNm = getWSprefix(inputWS) + 'conv_' + ftype + bgd + str(specMin) + "_to_" + str(specMax)
203    if Verbose:
204        logger.notice(func) 
205    PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func, 
206        StartX=startX, EndX=endX, FitType='Sequential')
207    wsname = confitParsToWS(outNm, inputWS, specMin, specMax)
208
209    # Add some information about convfit to the output workspace
210    options = getConvFitOption(ftype, bgd[:-2], Verbose)
211    AddSampleLog(Workspace=wsname, LogName="Fit Program", LogType="String", LogText='ConvFit')
212    AddSampleLog(Workspace=wsname, LogName='Background', LogType='String', LogText=str(options[0]))
213    AddSampleLog(Workspace=wsname, LogName='Delta', LogType='String', LogText=str(options[1]))
214    AddSampleLog(Workspace=wsname, LogName='Lorentzians', LogType='String', LogText=str(options[2]))
215
216    RenameWorkspace(InputWorkspace=outNm, OutputWorkspace=outNm + "_Parameters")
217    getConvFitResult(inputWS, resFile, outNm, ftype, bgd, specMin, specMax, Verbose)
218    if Save:
219        o_path = os.path.join(workdir, wsname+'.nxs')                    # path name for nxs file
220        if Verbose:
221            logger.notice('Creating file : '+o_path)
222        SaveNexusProcessed(InputWorkspace=wsname, Filename=o_path)
223    if Plot != 'None':
224        confitPlotSeq(wsname, Plot)
225    EndTime('ConvFit')
226
227##############################################################################
228# Elwin
229##############################################################################
230
231def GetTemperature(root,tempWS,log_type,Verbose):
232    (instr, run) = getInstrRun(root)
233    run_name = instr+run
234    log_name = run_name+'_'+log_type
235    run = mtd[tempWS].getRun()
236    unit1 = 'Temperature'               # default values
237    unit2 = 'K'
238    if log_type in run:                 # test logs in WS
239        tmp = run[log_type].value
240        temp = tmp[len(tmp)-1]
241        xval = temp
242        mess = ' Run : '+run_name +' ; Temperature in log = '+str(temp)
243    else:                               # logs not in WS
244        logger.notice('Log parameter not found')
245        log_file = log_name+'.txt'
246        log_path = FileFinder.getFullPath(log_file)
247        if (log_path == ''):            # log file does not exists
248            mess = ' Run : '+run_name +' ; Temperature file not found'
249            xval = int(run_name[-3:])
250            unit1 = 'Run-number'
251            unit2 = 'last 3 digits'
252        else:                           # get from log file
253            LoadLog(Workspace=tempWS, Filename=log_path)
254            run_logs = mtd[tempWS].getRun()
255            tmp = run_logs[log_name].value
256            temp = tmp[len(tmp)-1]
257            xval = temp
258            mess = ' Run : '+run_name+' ; Temperature in file = '+str(temp)
259    if Verbose:
260        logger.notice(mess)
261    unit = [unit1,unit2]
262    return xval,unit
263
264def elwin(inputFiles, eRange, log_type='sample', Normalise = False,
265        Save=False, Verbose=False, Plot=False): 
266    StartTime('ElWin')
267    workdir = config['defaultsave.directory']
268    CheckXrange(eRange,'Energy')
269    tempWS = '__temp'
270    Range2 = ( len(eRange) == 4 )
271    if Verbose:
272        range1 = str(eRange[0])+' to '+str(eRange[1])
273        if ( len(eRange) == 4 ): 
274            range2 = str(eRange[2])+' to '+str(eRange[3])
275            logger.notice('Using 2 energy ranges from '+range1+' & '+range2)
276        elif ( len(eRange) == 2 ):
277            logger.notice('Using 1 energy range from '+range1)
278    nr = 0
279    inputRuns = sorted(inputFiles)
280    for file in inputRuns:
281        (direct, file_name) = os.path.split(file)
282        (root, ext) = os.path.splitext(file_name)
283        LoadNexus(Filename=file, OutputWorkspace=tempWS)
284        nsam,ntc = CheckHistZero(tempWS)
285        (xval, unit) = GetTemperature(root,tempWS,log_type,Verbose)
286        if Verbose:
287            logger.notice('Reading file : '+file)
288        if ( len(eRange) == 4 ):
289            ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1], 
290                Range2Start=eRange[2], Range2End=eRange[3],
291                OutputInQ='__eq1', OutputInQSquared='__eq2')
292        elif ( len(eRange) == 2 ):
293            ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1],
294                OutputInQ='__eq1', OutputInQSquared='__eq2')
295        (instr, last) = getInstrRun(root)
296        q1 = np.array(mtd['__eq1'].readX(0))
297        i1 = np.array(mtd['__eq1'].readY(0))
298        e1 = np.array(mtd['__eq1'].readE(0))
299        Logarithm(InputWorkspace='__eq2', OutputWorkspace='__eq2')
300        q2 = np.array(mtd['__eq2'].readX(0))
301        i2 = np.array(mtd['__eq2'].readY(0))
302        e2 = np.array(mtd['__eq2'].readE(0))
303        if (nr == 0):
304            CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elf')
305            first = getWSprefix(tempWS,root)
306            datX1 = q1
307            datY1 = i1
308            datE1 = e1
309            datX2 = q2
310            datY2 = i2
311            datE2 = e2
312            Tvalue = [xval]
313            Terror = [0.0]
314            Taxis = str(xval)
315        else:
316            CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elftmp')
317            ConjoinWorkspaces(InputWorkspace1='__elf', InputWorkspace2='__elftmp', CheckOverlapping=False)
318            datX1 = np.append(datX1,q1)
319            datY1 = np.append(datY1,i1)
320            datE1 = np.append(datE1,e1)
321            datX2 = np.append(datX2,q2)
322            datY2 = np.append(datY2,i2)
323            datE2 = np.append(datE2,e2)
324            Tvalue.append(xval)
325            Terror.append(0.0)
326            Taxis += ','+str(xval)
327        nr += 1
328    Txa = np.array(Tvalue)
329    Tea = np.array(Terror)
330    nQ = len(q1)
331    for nq in range(0,nQ):
332        iq = []
333        eq = []
334        for nt in range(0,len(Tvalue)):
335            ii = mtd['__elf'].readY(nt)
336            iq.append(ii[nq])
337            ie = mtd['__elf'].readE(nt)
338            eq.append(ie[nq])
339        iqa = np.array(iq)
340        eqa = np.array(eq)
341        if (nq == 0):
342            datTx = Txa
343            datTy = iqa
344            datTe = eqa
345        else:
346            datTx = np.append(datTx,Txa)
347            datTy = np.append(datTy,iqa)
348            datTe = np.append(datTe,eqa)
349    DeleteWorkspace(tempWS)
350    DeleteWorkspace('__eq1')
351    DeleteWorkspace('__eq2')
352    if (nr == 1):
353        ename = first[:-1]
354    else:
355        ename = first+'to_'+last
356   
357    elfWS = ename+'_elf'    # interchange Q & T
358
359    #check if temp was increasing of decreasing
360    if(datTx[0] > datTx[-1]):
361        # if so reverse data to follow natural ordering
362        datTx = datTx[::-1]
363        datTy = datTy[::-1]
364        datTe = datTe[::-1]
365   
366    CreateWorkspace(OutputWorkspace=elfWS, DataX=datTx, DataY=datTy, DataE=datTe,
367        Nspec=nQ, UnitX='Energy', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q1)
368    DeleteWorkspace('__elf')
369    label = unit[0]+' / '+unit[1]
370    AddSampleLog(Workspace=elfWS, LogName="Vaxis", LogType="String", LogText=label)
371    e1WS = ename+'_eq1'
372    CreateWorkspace(OutputWorkspace=e1WS, DataX=datX1, DataY=datY1, DataE=datE1,
373        Nspec=nr, UnitX='MomentumTransfer', VerticalAxisUnit='Energy', VerticalAxisValues=Taxis)
374    label = unit[0]+' / '+unit[1]
375    AddSampleLog(Workspace=e1WS, LogName="Vaxis", LogType="String", LogText=label)
376    e2WS = ename+'_eq2'
377    CreateWorkspace(OutputWorkspace=e2WS, DataX=datX2, DataY=datY2, DataE=datE2,
378        Nspec=nr, UnitX='QSquared', VerticalAxisUnit='Energy', VerticalAxisValues=Taxis)
379    AddSampleLog(Workspace=e2WS, LogName="Vaxis", LogType="String", LogText=label)
380
381    if unit[0] == 'Temperature':
382        nT = len(Tvalue)
383        if Tvalue[0] < Tvalue[nT-1]:
384            lo = 0
385            hi = nT-1
386        else:
387            lo = nT-1
388            hi = 0
389        text = 'Temperature range : '+str(Tvalue[lo])+' to '+str(Tvalue[hi])
390        AddSampleLog(Workspace=e1WS, LogName="Temperature normalise", LogType="String", LogText=str(Normalise))
391        if Normalise:
392            yval = mtd[e1WS].readY(lo)
393            normFactor = 1.0/yval[0]
394            Scale(InputWorkspace=e1WS, OutputWorkspace=e1WS, Factor=normFactor, Operation='Multiply')
395            AddSampleLog(Workspace=e1WS, LogName="Temperature value", LogType="Number", LogText=str(Tvalue[0]))
396            if Verbose:
397                text = 'Temperature range : '+str(Tvalue[lo])+' to '+str(Tvalue[hi])
398                logger.notice(text)
399                logger.notice('Normalised eq1 by scale factor : '+str(normFactor))
400
401    unity = mtd[e1WS].getAxis(1).setUnit("Label")
402    unity.setLabel(unit[0], unit[1])
403    label = unit[0]+' / '+unit[1]
404    addElwinLogs(e1WS, label, eRange, Range2)
405   
406    unity = mtd[e2WS].getAxis(1).setUnit("Label")
407    unity.setLabel(unit[0], unit[1])
408    addElwinLogs(e2WS, label, eRange, Range2)
409   
410    unitx = mtd[elfWS].getAxis(0).setUnit("Label")
411    unitx.setLabel(unit[0], unit[1])
412    addElwinLogs(elfWS, label, eRange, Range2)
413
414    if Save:
415        e1_path = os.path.join(workdir, e1WS+'.nxs')                                    # path name for nxs file
416        e2_path = os.path.join(workdir, e2WS+'.nxs')                                    # path name for nxs file
417        elf_path = os.path.join(workdir, elfWS+'.nxs')                                  # path name for nxs file
418
419        if Verbose:
420            logger.notice('Creating file : '+e1_path)
421            logger.notice('Creating file : '+e2_path)
422            logger.notice('Creating file : '+elf_path)
423
424        SaveNexusProcessed(InputWorkspace=e1WS, Filename=e1_path)
425        SaveNexusProcessed(InputWorkspace=e2WS, Filename=e2_path)
426        SaveNexusProcessed(InputWorkspace=elfWS, Filename=elf_path)
427
428    if Plot:
429        elwinPlot(label,e1WS,e2WS,elfWS)
430
431    EndTime('Elwin')
432    return e1WS,e2WS
433
434# Add sample log to each of the workspaces created by Elwin
435def addElwinLogs(ws, label, eRange, Range2):
436
437    AddSampleLog(Workspace=ws, LogName="Vaxis", LogType="String", LogText=label)
438    AddSampleLog(Workspace=ws, LogName="Range1 start", LogType="Number", LogText=str(eRange[0]))
439    AddSampleLog(Workspace=ws, LogName="Range1 end", LogType="Number", LogText=str(eRange[1]))
440    AddSampleLog(Workspace=ws, LogName="Two ranges", LogType="String", LogText=str(Range2))
441
442    if Range2:
443        AddSampleLog(Workspace=ws, LogName="Range2 start", LogType="Number", LogText=str(eRange[2]))
444        AddSampleLog(Workspace=ws, LogName="Range2 end", LogType="Number", LogText=str(eRange[3]))
445
446#Plot each of the workspace output by elwin
447def elwinPlot(label,eq1,eq2,elf):
448    plotElwinWorkspace(eq1, yAxisTitle='Elastic Intensity', setScale=True)
449    plotElwinWorkspace(eq2, yAxisTitle='log(Elastic Intensity)', setScale=True)
450    plotElwinWorkspace(elf, xAxisTitle=label)
451
452#Plot a workspace generated by Elwin
453def plotElwinWorkspace(ws, xAxisTitle=None, yAxisTitle=None, setScale=False):
454    ws = mtd[ws]
455    nBins = ws.blocksize()
456    lastX = ws.readX(0)[nBins-1]
457
458    nhist = ws.getNumberHistograms()
459
460    try:
461        graph = mp.plotSpectrum(ws, range(0,nhist))
462    except RuntimeError, e:
463        #User clicked cancel on plot so don't do anything
464        return None
465   
466    layer = graph.activeLayer()
467
468    #set the x scale of the layer
469    if setScale:
470        layer.setScale(mp.Layer.Bottom, 0.0, lastX)
471   
472    #set the title on the x-axis
473    if xAxisTitle:
474        layer.setAxisTitle(mp.Layer.Bottom, xAxisTitle)
475
476    #set the title on the y-axis
477    if yAxisTitle:
478        layer.setAxisTitle(mp.Layer.Left, yAxisTitle)
479
480##############################################################################
481# Fury
482##############################################################################
483
484def furyPlot(inWS, spec):
485    graph = mp.plotSpectrum(inWS, spec)
486    layer = graph.activeLayer()
487    layer.setScale(mp.Layer.Left, 0, 1.0)
488
489def fury(samWorkspaces, res_file, rebinParam, RES=True, Save=False, Verbose=False,
490        Plot=False): 
491    StartTime('Fury')
492    workdir = config['defaultsave.directory']
493    samTemp = samWorkspaces[0]
494    nsam,npt = CheckHistZero(samTemp)
495    Xin = mtd[samTemp].readX(0)
496    d1 = Xin[1]-Xin[0]
497    if d1 < 1e-8:
498        error = 'Data energy bin is zero'
499        logger.notice('ERROR *** ' + error)
500        sys.exit(error)
501    d2 = Xin[npt-1]-Xin[npt-2]
502    dmin = min(d1,d2)
503    pars = rebinParam.split(',')
504    if (float(pars[1]) <= dmin):
505        error = 'EWidth = ' + pars[1] + ' < smallest Eincr = ' + str(dmin)
506        logger.notice('ERROR *** ' + error)
507        sys.exit(error)
508    outWSlist = []
509    # Process RES Data Only Once
510    if Verbose:
511        logger.notice('Reading RES file : '+res_file)
512    LoadNexus(Filename=res_file, OutputWorkspace='res_data') # RES
513    CheckAnalysers(samTemp,'res_data',Verbose)
514    nres,nptr = CheckHistZero('res_data')
515    if nres > 1:
516        CheckHistSame(samTemp,'Sample','res_data','Resolution')
517    Rebin(InputWorkspace='res_data', OutputWorkspace='res_data', Params=rebinParam)
518    Integration(InputWorkspace='res_data', OutputWorkspace='res_int')
519    ConvertToPointData(InputWorkspace='res_data', OutputWorkspace='res_data')
520    ExtractFFTSpectrum(InputWorkspace='res_data', OutputWorkspace='res_fft', FFTPart=2)
521    Divide(LHSWorkspace='res_fft', RHSWorkspace='res_int', OutputWorkspace='res')
522    for samWs in samWorkspaces:
523        (direct, filename) = os.path.split(samWs)
524        (root, ext) = os.path.splitext(filename)
525        Rebin(InputWorkspace=samWs, OutputWorkspace='sam_data', Params=rebinParam)
526        Integration(InputWorkspace='sam_data', OutputWorkspace='sam_int')
527        ConvertToPointData(InputWorkspace='sam_data', OutputWorkspace='sam_data')
528        ExtractFFTSpectrum(InputWorkspace='sam_data', OutputWorkspace='sam_fft', FFTPart=2)
529        Divide(LHSWorkspace='sam_fft', RHSWorkspace='sam_int', OutputWorkspace='sam')
530        # Create save file name
531        savefile = getWSprefix('sam_data', root) + 'iqt'
532        outWSlist.append(savefile)
533        Divide(LHSWorkspace='sam', RHSWorkspace='res', OutputWorkspace=savefile)
534        #Cleanup Sample Files
535        DeleteWorkspace('sam_data')
536        DeleteWorkspace('sam_int')
537        DeleteWorkspace('sam_fft')
538        DeleteWorkspace('sam')
539        # Crop nonsense values off workspace
540        bin = int(math.ceil(mtd[savefile].blocksize()/2.0))
541        binV = mtd[savefile].dataX(0)[bin]
542        CropWorkspace(InputWorkspace=savefile, OutputWorkspace=savefile, XMax=binV)
543        if Save:
544            opath = os.path.join(workdir, savefile+'.nxs')                                      # path name for nxs file
545            SaveNexusProcessed(InputWorkspace=savefile, Filename=opath)
546            if Verbose:
547                logger.notice('Output file : '+opath) 
548    # Clean Up RES files
549    DeleteWorkspace('res_data')
550    DeleteWorkspace('res_int')
551    DeleteWorkspace('res_fft')
552    DeleteWorkspace('res')
553    if Plot:
554        specrange = range(0,mtd[outWSlist[0]].getNumberHistograms())
555        furyPlot(outWSlist, specrange)
556    EndTime('Fury')
557    return outWSlist
558
559##############################################################################
560# FuryFit
561##############################################################################
562
563def getFuryFitOption(option):
564    nopt = len(option)
565    if nopt == 2:
566        npeak = option[0]
567        type = option[1]
568    elif nopt == 4:
569        npeak = '2'
570        type = 'SE'
571    else:
572        error = 'Bad option : ' +option
573        logger.notice('ERROR *** ' + error)
574        sys.exit(error)
575    return npeak, type
576
577def furyfitParsToWS(Table, Data, option):
578    npeak, type = getFuryFitOption(option)   
579    Q = createQaxis(Data)
580    nQ = len(Q)
581    ws = mtd[Table]
582    rCount = ws.rowCount()
583    cCount = ws.columnCount()
584    cName =  ws.getColumnNames()
585    Qa = np.array(Q)
586    A0v = ws.column(1)     #bgd value
587    A0e = ws.column(2)     #bgd error
588    Iy1 = ws.column(5)      #intensity1 value
589    Ie1 = ws.column(2)      #intensity1 error = bgd
590    dataX = Qa
591    dataY = np.array(A0v)
592    dataE = np.array(A0e)
593    names = cName[1]
594    dataX = np.append(dataX,Qa)
595    dataY = np.append(dataY,np.array(Iy1))
596    dataE = np.append(dataE,np.array(Ie1))
597    names += ","+cName[5]
598    Ty1 = ws.column(7)      #tau1 value
599    Te1 = ws.column(8)      #tau1 error
600    dataX = np.append(dataX,Qa)
601    dataY = np.append(dataY,np.array(Ty1))
602    dataE = np.append(dataE,np.array(Te1))
603    names += ","+cName[7]
604    nSpec = 3
605    if npeak == '1' and type == 'S':
606        By1 = ws.column(9)  #beta1 value
607        Be1 = ws.column(10) #beta2 error
608        dataX = np.append(dataX,Qa)
609        dataY = np.append(dataY,np.array(By1))
610        dataE = np.append(dataE,np.array(Be1))
611        names += ","+cName[9]
612        nSpec += 1
613    if npeak == '2':
614        Iy2 = ws.column(9)  #intensity2 value
615        Ie2 = ws.column(10) #intensity2 error
616        dataX = np.append(dataX,Qa)
617        dataY = np.append(dataY,np.array(Iy2))
618        dataE = np.append(dataE,np.array(Ie2))
619        names += ","+cName[9]
620        nSpec += 1
621        Ty2 = ws.column(11)  #tau2 value
622        Te2 = ws.column(12) #tau2 error
623        dataX = np.append(dataX,Qa)
624        dataY = np.append(dataY,np.array(Ty2))
625        dataE = np.append(dataE,np.array(Te2))
626        names += ","+cName[11]
627        nSpec += 1
628    wsname = Table + "_Workspace"
629    CreateWorkspace(OutputWorkspace=wsname, DataX=dataX, DataY=dataY, DataE=dataE, 
630        Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
631        VerticalAxisValues=names)
632    return wsname
633
634def createFurySeqResFun(ties, par, option):
635    npeak, type = getFuryFitOption(option)   
636    fun = 'name=LinearBackground,A0='+str(par[0])+',A1=0,ties=(A1=0);'
637    if npeak == '1':
638        if type == 'E':
639            fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)),Intensity='+str(par[1])+',Tau='+str(par[2])
640        if type == 'S':
641            fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity='+str(par[1])+',Tau='+str(par[2])+',Beta='+str(par[3])
642        if ties:
643            fun += ';ties=(f1.Intensity=1-f0.A0)'
644    if npeak == '2':
645        fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)),Intensity='+str(par[1])+',Tau='+str(par[2])
646        if type == 'E':
647            fun += ';name=UserFunction,Formula=Intensity*exp(-(x/Tau)),Intensity='+str(par[3])+',Tau='+str(par[4])
648        if type == 'SE':
649            fun += ';name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity='+str(par[3])+',Tau='+str(par[4])+',Beta='+str(par[5])
650        if ties:
651            fun += ';ties=(f1.Intensity=1-f2.Intensity-f0.A0)'
652    return fun
653
654def getFurySeqResult(inputWS, outNm, option, Verbose):
655    logger.notice('Option : ' +option)
656    npeak, type = getFuryFitOption(option)   
657    params = mtd[outNm+'_Parameters']
658    A0 = params.column(1)     #bgd value
659    I1 = params.column(5)      #intensity1 value
660    T1 = params.column(7)      #tau1 value
661    if npeak == '1' and type == 'S':
662        B1 = params.column(9)  #beta1 value
663    if npeak == '2':
664        I2 = params.column(9)  #intensity2 value
665        T2 = params.column(11)  #tau2 value
666        if type == 'SE':
667            B2 = params.column(13)  #beta2 value
668    nHist = mtd[inputWS].getNumberHistograms()
669    for i in range(nHist):
670        paras = [A0[i], I1[i], T1[i]]
671        if npeak == '1' and type == 'S':
672            paras.append(B1[i])
673        if npeak == '2':
674            paras.append(I2[i])
675            paras.append(T2[i])
676            if type == 'SE':
677                paras.append(B2[i])
678        func = createFurySeqResFun(True, paras, option)
679        if Verbose:
680            logger.notice('Fit func : '+func)   
681        fitWS = outNm + '_Result_'
682        fout = fitWS + str(i)
683        Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i,Output=fout,MaxIterations=0)
684        unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
685        unitx.setLabel('Time' , 'ns')
686        RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
687        DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
688        DeleteWorkspace(fitWS+str(i)+'_Parameters')
689        if i == 0:
690            group = fout
691        else:
692            group += ',' + fout
693    GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
694
695def furyfitPlotSeq(inputWS, Plot):
696    nHist = mtd[inputWS].getNumberHistograms()
697    if ( Plot == 'All' ):
698        mp.plotSpectrum(inputWS, range(0, nHist), True)
699        return
700    plotSpecs = []
701    if ( Plot == 'Intensity' ):
702        res = 'Intensity$'
703    if ( Plot == 'Tau' ):
704        res = 'Tau$'
705    elif ( Plot == 'Beta' ):
706        res = 'Beta$'   
707    for i in range(0, nHist):
708        title = mtd[inputWS].getAxis(1).label(i)
709        if ( re.search(res, title) ):
710            plotSpecs.append(i)
711    mp.plotSpectrum(inputWS, plotSpecs, True)
712
713def furyfitSeq(inputWS, func, ftype, startx, endx, Save, Plot, Verbose=False): 
714    StartTime('FuryFit')
715    workdir = config['defaultsave.directory']
716    input = inputWS+',i0'
717    nHist = mtd[inputWS].getNumberHistograms()
718    for i in range(1,nHist):
719        input += ';'+inputWS+',i'+str(i)
720    outNm = getWSprefix(inputWS) + 'fury_' + ftype + "0_to_" + str(nHist-1)
721    option = ftype[:-2]
722    if Verbose:
723        logger.notice('Option: '+option) 
724        logger.notice(func) 
725    PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func, 
726        StartX=startx, EndX=endx, FitType='Sequential')
727    fitWS = furyfitParsToWS(outNm, inputWS, option)
728    RenameWorkspace(InputWorkspace=outNm, OutputWorkspace=outNm+"_Parameters")
729    CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=inputWS, XMin=startx, XMax=endx)
730    getFurySeqResult(inputWS, outNm, option, Verbose)
731    CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fitWS)
732    AddSampleLog(Workspace=fitWS, LogName="start_x", LogType="Number", LogText=str(startx))
733    AddSampleLog(Workspace=fitWS, LogName="end_x", LogType="Number", LogText=str(endx))
734    AddSampleLog(Workspace=fitWS, LogName="option", LogType="String", LogText=option)
735    CopyLogs(InputWorkspace=inputWS, OutputWorkspace=outNm+"_Result")
736    AddSampleLog(Workspace=outNm+'_Result', LogName="start_x", LogType="Number", LogText=str(startx))
737    AddSampleLog(Workspace=outNm+'_Result', LogName="end_x", LogType="Number", LogText=str(endx))
738    AddSampleLog(Workspace=outNm+'_Result', LogName="option", LogType="String", LogText=option)
739    if Save:
740        fpath = os.path.join(workdir, fitWS+'.nxs')                                     # path name for nxs file
741        SaveNexusProcessed(InputWorkspace=fitWS, Filename=fpath)
742        rpath = os.path.join(workdir, outNm+'_Result.nxs')                                      # path name for nxs file
743        SaveNexusProcessed(InputWorkspace=outNm+'_Result', Filename=fpath)
744        if Verbose:
745            logger.notice('Fit output file : '+fpath) 
746            logger.notice('Result output file : '+rpath) 
747    if ( Plot != 'None' ):
748        furyfitPlotSeq(fitWS, Plot)
749    EndTime('FuryFit')
750    return mtd[fitWS]
751
752def furyfitMultParsToWS(Table, Data):
753#   Q = createQaxis(Data)
754    theta,Q = GetThetaQ(Data)
755    ws = mtd[Table+'_Parameters']
756    rCount = ws.rowCount()
757    cCount = ws.columnCount()
758    nSpec = ( rCount - 1 ) / 5
759    val = ws.column(1)     #value
760    err = ws.column(2)     #error
761    dataX = []
762    A0val = []
763    A0err = []
764    Ival = []
765    Ierr = []
766    Tval = []
767    Terr = []
768    Bval = []
769    Berr = []
770    for spec in range(0,nSpec):
771        n1 = spec*5
772        A0 = n1
773        A1 = n1+1
774        int = n1+2                   #intensity value
775        tau = n1+3                   #tau value
776        beta = n1+4                   #beta value
777        dataX.append(Q[spec])
778        A0val.append(val[A0])
779        A0err.append(err[A0])
780        Ival.append(val[int])
781        Ierr.append(err[int])
782        Tval.append(val[tau])
783        Terr.append(err[tau])
784        Bval.append(val[beta])
785        Berr.append(err[beta])
786    nQ = len(dataX)
787    Qa = np.array(dataX)
788    dataY = np.array(A0val)
789    dataE = np.array(A0err)
790    dataY = np.append(dataY,np.array(Ival))
791    dataE = np.append(dataE,np.array(Ierr))
792    dataY = np.append(dataY,np.array(Tval))
793    dataE = np.append(dataE,np.array(Terr))
794    dataY = np.append(dataY,np.array(Bval))
795    dataE = np.append(dataE,np.array(Berr))
796    names = 'A0,Intensity,Tau,Beta'
797    suffix = 'Workspace'
798    wsname = Table + '_' + suffix
799    CreateWorkspace(OutputWorkspace=wsname, DataX=Qa, DataY=dataY, DataE=dataE, 
800        Nspec=4, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
801        VerticalAxisValues=names)
802    return wsname
803
804def furyfitPlotMult(inputWS, Plot):
805    nHist = mtd[inputWS].getNumberHistograms()
806    if ( Plot == 'All' ):
807        mp.plotSpectrum(inputWS, range(0, nHist))
808        return
809    plotSpecs = []
810    if ( Plot == 'Intensity' ):
811        mp.plotSpectrum(inputWS, 1, True)
812    if ( Plot == 'Tau' ):
813        mp.plotSpectrum(inputWS, 2, True)
814    elif ( Plot == 'Beta' ):
815        mp.plotSpectrum(inputWS, 3, True)   
816
817
818def createFuryMultFun(ties = True, function = ''):
819    fun =  '(composite=CompositeFunction,$domains=i;'
820    fun += function
821    if ties:
822        fun += ';ties=(f1.Intensity=1-f0.A0)'
823    fun += ');'
824    return fun
825
826def createFuryMultResFun(ties = True, A0 = 0.02, Intensity = 0.98 ,Tau = 0.025, Beta = 0.8):
827    fun =  '(composite=CompositeFunction,$domains=i;'
828    fun += 'name=LinearBackground,A0='+str(A0)+',A1=0,ties=(A1=0);'
829    fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity='+str(Intensity)+',Tau='+str(Tau)+',Beta='+str(Beta)
830    if ties:
831        fun += ';ties=(f1.Intensity=1-f0.A0)'
832    fun += ');'
833    return fun
834
835def getFuryMultResult(inputWS, outNm, function, Verbose):
836    params = mtd[outNm+'_Parameters']
837    nHist = mtd[inputWS].getNumberHistograms()
838    for i in range(nHist):
839        j = 5 * i
840#        assert( params.row(j)['Name'][3:] == 'f0.A0' )
841        A0 = params.row(j)['Value']
842        A1 = params.row(j + 1)['Value']
843        Intensity = params.row(j + 2)['Value']
844        Tau = params.row(j + 3)['Value']
845        Beta = params.row(j + 4)['Value']
846        func = createFuryMultResFun(True,  A0, Intensity ,Tau, Beta)
847        if Verbose:
848            logger.notice('Fit func : '+func)   
849        fitWS = outNm + '_Result_'
850        fout = fitWS + str(i)
851        Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i,Output=fout,MaxIterations=0)
852        unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
853        unitx.setLabel('Time' , 'ns')
854        RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
855        DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
856        DeleteWorkspace(fitWS+str(i)+'_Parameters')
857        if i == 0:
858            group = fout
859        else:
860            group += ',' + fout
861    GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
862
863def furyfitMult(inputWS, function, ftype, startx, endx, Save, Plot, Verbose=False):
864    StartTime('FuryFit Mult')
865    workdir = config['defaultsave.directory']
866    option = ftype[:-2]
867    if Verbose:
868        logger.notice('Option: '+option) 
869        logger.notice('Function: '+function) 
870    nHist = mtd[inputWS].getNumberHistograms()
871    outNm = inputWS[:-3] + 'fury_mult'
872    f1 = createFuryMultFun(True, function)
873    func= 'composite=MultiDomainFunction,NumDeriv=1;'
874    ties='ties=('
875    kwargs = {}
876    for i in range(0,nHist):
877        func+=f1
878        if i > 0:
879            ties += 'f' + str(i) + '.f1.Beta=f0.f1.Beta'
880            if i < nHist-1:
881                ties += ','
882            kwargs['InputWorkspace_' + str(i)] = inputWS
883        kwargs['WorkspaceIndex_' + str(i)] = i
884    ties+=')'
885    func += ties
886    CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=inputWS, XMin=startx, XMax=endx)
887    Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=0,Output=outNm,**kwargs)
888    outWS = furyfitMultParsToWS(outNm, inputWS)
889    getFuryMultResult(inputWS, outNm, function, Verbose)
890    CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fitWS)
891    AddSampleLog(Workspace=fitWS, LogName="start_x", LogType="Number", LogText=str(startx))
892    AddSampleLog(Workspace=fitWS, LogName="end_x", LogType="Number", LogText=str(endx))
893    AddSampleLog(Workspace=fitWS, LogName="function", LogType="String", LogText=ftype)
894    AddSampleLog(Workspace=fitWS, LogName="link_q", LogType="String", LogText='True')
895    CopyLogs(InputWorkspace=inputWS, OutputWorkspace=outNm+'_result')
896    AddSampleLog(Workspace=outNm+'_result', LogName="start_x", LogType="Number", LogText=str(startx))
897    AddSampleLog(Workspace=outNm+'_result', LogName="end_x", LogType="Number", LogText=str(endx))
898    AddSampleLog(Workspace=outNm+'_result', LogName="function", LogType="String", LogText=ftype)
899    AddSampleLog(Workspace=outNm+'_result', LogName="link_q", LogType="String", LogText='True')
900    if Save:
901        opath = os.path.join(workdir, outWS+'.nxs')                                     # path name for nxs file
902        SaveNexusProcessed(InputWorkspace=outWS, Filename=opath)
903        rpath = os.path.join(workdir, outNm+'_result.nxs')                                      # path name for nxs file
904        SaveNexusProcessed(InputWorkspace=outNm+'_result', Filename=rpath)
905        if Verbose:
906            logger.notice('Output file : '+opath) 
907            logger.notice('Output file : '+rpath) 
908    if ( Plot != 'None' ):
909        furyfitPlotMult(outWS, Plot)
910    EndTime('FuryFit')
911
912##############################################################################
913# MSDFit
914##############################################################################
915
916def msdfitParsToWS(Table, xData):
917    dataX = xData
918    ws = mtd[Table+'_Table']
919    rCount = ws.rowCount()
920    yA0 = ws.column(1)
921    eA0 = ws.column(2)
922    yA1 = ws.column(3) 
923    dataY1 = map(lambda x : -x, yA1) 
924    eA1 = ws.column(4)
925    wsname = Table
926    CreateWorkspace(OutputWorkspace=wsname+'_a0', DataX=dataX, DataY=yA0, DataE=eA0,
927        Nspec=1, UnitX='')
928    CreateWorkspace(OutputWorkspace=wsname+'_a1', DataX=dataX, DataY=dataY1, DataE=eA1,
929        Nspec=1, UnitX='')
930    group = wsname+'_a0,'+wsname+'_a1'
931    GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=wsname)
932    return wsname
933
934def msdfitPlotSeq(inputWS, xlabel):
935    msd_plot = mp.plotSpectrum(inputWS+'_a1',0,True)
936    msd_layer = msd_plot.activeLayer()
937    msd_layer.setAxisTitle(mp.Layer.Bottom,xlabel)
938    msd_layer.setAxisTitle(mp.Layer.Left,'<u2>')
939
940def msdfitPlotFits(calcWS, n):
941    mfit_plot = mp.plotSpectrum(calcWS+'_0',[0,1],True)
942    mfit_layer = mfit_plot.activeLayer()
943    mfit_layer.setAxisTitle(mp.Layer.Left,'log(Elastic Intensity)')
944
945def msdfit(inputs, startX, endX, Save=False, Verbose=False, Plot=True): 
946    StartTime('msdFit')
947    workdir = config['defaultsave.directory']
948    log_type = 'sample'
949    file = inputs[0]
950    (direct, filename) = os.path.split(file)
951    (root, ext) = os.path.splitext(filename)
952    (instr, first) = getInstrRun(filename)
953    if Verbose:
954        logger.notice('Reading Run : '+file)
955    LoadNexusProcessed(FileName=file, OutputWorkspace=root)
956    nHist = mtd[root].getNumberHistograms()
957    file_list = []
958    run_list = []
959    ws = mtd[root]
960    ws_run = ws.getRun()
961    vertAxisValues = ws.getAxis(1).extractValues()
962    x_list = vertAxisValues
963    if 'Vaxis' in ws_run:
964        xlabel = ws_run.getLogData('Vaxis').value
965    for nr in range(0, nHist):
966        nsam,ntc = CheckHistZero(root)
967        lnWS = '__lnI_'+str(nr)
968        file_list.append(lnWS)
969        ExtractSingleSpectrum(InputWorkspace=root, OutputWorkspace=lnWS,
970            WorkspaceIndex=nr)
971        if (nr == 0):
972            run_list = lnWS
973        else:
974            run_list += ';'+lnWS
975    mname = root[:-4]
976    msdWS = mname+'_msd'
977    if Verbose:
978       logger.notice('Fitting Runs '+mname)
979       logger.notice('Q-range from '+str(startX)+' to '+str(endX))
980    function = 'name=LinearBackground, A0=0, A1=0'
981    PlotPeakByLogValue(Input=run_list, OutputWorkspace=msdWS+'_Table', Function=function,
982        StartX=startX, EndX=endX, FitType = 'Sequential')
983    msdfitParsToWS(msdWS, x_list)
984    nr = 0
985    fitWS = mname+'_Fit'
986    calcWS = mname+'_msd_Result'
987    a0 = mtd[msdWS+'_a0'].readY(0)
988    a1 = mtd[msdWS+'_a1'].readY(0)
989    for nr in range(0, nHist):
990        inWS = file_list[nr]
991        CropWorkspace(InputWorkspace=inWS,OutputWorkspace='__data',XMin=0.95*startX,XMax=1.05*endX)
992        dataX = mtd['__data'].readX(0)
993        nxd = len(dataX)
994        dataX = np.append(dataX,2*dataX[nxd-1]-dataX[nxd-2])
995        dataY = np.array(mtd['__data'].readY(0))
996        dataE = np.array(mtd['__data'].readE(0))
997        xd = []
998        yd = []
999        ed = []
1000        for n in range(0,nxd):
1001            line = a0[nr] - a1[nr]*dataX[n]
1002            xd.append(dataX[n])
1003            yd.append(line)
1004            ed.append(0.0)
1005        xd.append(dataX[nxd])
1006        dataX = np.append(dataX,np.array(xd))
1007        dataY = np.append(dataY,np.array(yd))
1008        dataE = np.append(dataE,np.array(ed))
1009        fout = calcWS +'_'+ str(nr)
1010        CreateWorkspace(OutputWorkspace=fout, DataX=dataX, DataY=dataY, DataE=dataE,
1011            Nspec=2, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues='Data,Calc')
1012        if nr == 0:
1013            gro = fout
1014        else:
1015            gro += ',' + fout
1016        DeleteWorkspace(inWS)
1017        DeleteWorkspace('__data')
1018    GroupWorkspaces(InputWorkspaces=gro,OutputWorkspace=calcWS)
1019    if Plot:
1020        msdfitPlotSeq(msdWS, xlabel)
1021        msdfitPlotFits(calcWS, 0)
1022    if Save:
1023        msd_path = os.path.join(workdir, msdWS+'.nxs')                                  # path name for nxs file
1024        SaveNexusProcessed(InputWorkspace=msdWS, Filename=msd_path, Title=msdWS)
1025        if Verbose:
1026            logger.notice('Output msd file : '+msd_path) 
1027    EndTime('msdFit')
1028    return msdWS
1029
1030def plotInput(inputfiles,spectra=[]):
1031    OneSpectra = False
1032    if len(spectra) != 2:
1033        spectra = [spectra[0], spectra[0]]
1034        OneSpectra = True
1035    workspaces = []
1036    for file in inputfiles:
1037        root = LoadNexus(Filename=file)
1038        if not OneSpectra:
1039            GroupDetectors(root, root,
1040                DetectorList=range(spectra[0],spectra[1]+1) )
1041        workspaces.append(root)
1042    if len(workspaces) > 0:
1043        graph = mp.plotSpectrum(workspaces,0)
1044        layer = graph.activeLayer().setTitle(", ".join(workspaces))
1045       
1046##############################################################################
1047# Corrections
1048##############################################################################
1049
1050def CubicFit(inputWS, spec, Verbose=False):
1051    '''Uses the Mantid Fit Algorithm to fit a quadratic to the inputWS
1052    parameter. Returns a list containing the fitted parameter values.'''
1053    function = 'name=Quadratic, A0=1, A1=0, A2=0'
1054    fit = Fit(Function=function, InputWorkspace=inputWS, WorkspaceIndex=spec,
1055      CreateOutput=True, Output='Fit')
1056    table = mtd['Fit_Parameters']
1057    A0 = table.cell(0,1)
1058    A1 = table.cell(1,1)
1059    A2 = table.cell(2,1)
1060    Abs = [A0, A1, A2]
1061    if Verbose:
1062       logger.notice('Group '+str(spec)+' of '+inputWS+' ; fit coefficients are : '+str(Abs))
1063    return Abs
1064
1065def applyCorrections(inputWS, canWS, corr, Verbose=False):
1066    '''Through the PolynomialCorrection algorithm, makes corrections to the
1067    input workspace based on the supplied correction values.'''
1068    # Corrections are applied in Lambda (Wavelength)
1069    efixed = getEfixed(inputWS)                # Get efixed
1070    theta,Q = GetThetaQ(inputWS)
1071    sam_name = getWSprefix(inputWS)
1072    ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='Wavelength',
1073        EMode='Indirect', EFixed=efixed)
1074
1075    nameStem = corr[:-4]
1076    if canWS != '':
1077        (instr, can_run) = getInstrRun(canWS)
1078        corrections = [nameStem+'_ass', nameStem+'_assc', nameStem+'_acsc', nameStem+'_acc']
1079        CorrectedWS = sam_name +'Correct_'+ can_run
1080        ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='Wavelength',
1081            EMode='Indirect', EFixed=efixed)
1082    else:
1083        corrections = [nameStem+'_ass']
1084        CorrectedWS = sam_name +'Corrected'
1085    nHist = mtd[inputWS].getNumberHistograms()
1086    # Check that number of histograms in each corrections workspace matches
1087    # that of the input (sample) workspace
1088    for ws in corrections:
1089        if ( mtd[ws].getNumberHistograms() != nHist ):
1090            raise ValueError('Mismatch: num of spectra in '+ws+' and inputWS')
1091    # Workspaces that hold intermediate results
1092    CorrectedSampleWS = '__csam'
1093    CorrectedCanWS = '__ccan'
1094    for i in range(0, nHist): # Loop through each spectra in the inputWS
1095        ExtractSingleSpectrum(InputWorkspace=inputWS, OutputWorkspace=CorrectedSampleWS,
1096            WorkspaceIndex=i)
1097        if ( len(corrections) == 1 ):
1098            Ass = CubicFit(corrections[0], i, Verbose)
1099            PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
1100                Coefficients=Ass, Operation='Divide')
1101            if ( i == 0 ):
1102                CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
1103            else:
1104                ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS)
1105        else:
1106            ExtractSingleSpectrum(InputWorkspace=canWS, OutputWorkspace=CorrectedCanWS,
1107                WorkspaceIndex=i)
1108            Acc = CubicFit(corrections[3], i, Verbose)
1109            PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
1110                Coefficients=Acc, Operation='Divide')
1111            Acsc = CubicFit(corrections[2], i, Verbose)
1112            PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
1113                Coefficients=Acsc, Operation='Multiply')
1114            Minus(LHSWorkspace=CorrectedSampleWS, RHSWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedSampleWS)
1115            Assc = CubicFit(corrections[1], i, Verbose)
1116            PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
1117                Coefficients=Assc, Operation='Divide')
1118            if ( i == 0 ):
1119                CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
1120            else:
1121                ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS,
1122                                  CheckOverlapping=False)
1123    ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='DeltaE',
1124        EMode='Indirect', EFixed=efixed)
1125    ConvertUnits(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS, Target='DeltaE',
1126        EMode='Indirect', EFixed=efixed)
1127    ConvertSpectrumAxis(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_rqw', 
1128        Target='ElasticQ', EMode='Indirect', EFixed=efixed)
1129    RenameWorkspace(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_red')
1130    if canWS != '':
1131        DeleteWorkspace(CorrectedCanWS)
1132        ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='DeltaE',
1133            EMode='Indirect', EFixed=efixed)
1134    DeleteWorkspace('Fit_NormalisedCovarianceMatrix')
1135    DeleteWorkspace('Fit_Parameters')
1136    DeleteWorkspace('Fit_Workspace')
1137    DeleteWorkspace(corr)
1138    return CorrectedWS
1139               
1140def abscorFeeder(sample, container, geom, useCor, corrections, Verbose=False, ScaleOrNotToScale=False, factor=1, Save=False,
1141        PlotResult='None', PlotContrib=False):
1142    '''Load up the necessary files and then passes them into the main
1143    applyCorrections routine.'''
1144    StartTime('ApplyCorrections')
1145    workdir = config['defaultsave.directory']
1146    s_hist,sxlen = CheckHistZero(sample)
1147    sam_name = getWSprefix(sample)
1148    efixed = getEfixed(sample)
1149    if container != '':
1150        CheckAnalysers(sample,container,Verbose)
1151        CheckHistSame(sample,'Sample',container,'Container')
1152        (instr, can_run) = getInstrRun(container)
1153        if ScaleOrNotToScale:
1154            Scale(InputWorkspace=container, OutputWorkspace=container, Factor=factor, Operation='Multiply')
1155            if Verbose:
1156                logger.notice('Container scaled by '+str(factor))
1157    if useCor:
1158        if Verbose:
1159            text = 'Correcting sample ' + sample
1160            if container != '':
1161                text += ' with ' + container
1162            logger.notice(text)
1163           
1164        cor_result = applyCorrections(sample, container, corrections, Verbose)
1165        rws = mtd[cor_result+'_red']
1166        outNm= cor_result + '_Result_'
1167
1168        if Save:
1169            cred_path = os.path.join(workdir,cor_result+'_red.nxs')
1170            SaveNexusProcessed(InputWorkspace=cor_result+'_red',Filename=cred_path)
1171            if Verbose:
1172                logger.notice('Output file created : '+cred_path)
1173        calc_plot = [cor_result+'_red',sample]
1174        res_plot = cor_result+'_rqw'
1175    else:
1176        if ( container == '' ):
1177            sys.exit('ERROR *** Invalid options - nothing to do!')
1178        else:
1179            sub_result = sam_name +'Subtract_'+ can_run
1180            if Verbose:
1181                logger.notice('Subtracting '+container+' from '+sample)
1182            Minus(LHSWorkspace=sample,RHSWorkspace=container,OutputWorkspace=sub_result)
1183            ConvertSpectrumAxis(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_rqw', 
1184                Target='ElasticQ', EMode='Indirect', EFixed=efixed)
1185            RenameWorkspace(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_red')
1186            rws = mtd[sub_result+'_red']
1187            outNm= sub_result + '_Result_'
1188            if Save:
1189                sred_path = os.path.join(workdir,sub_result+'_red.nxs')
1190                SaveNexusProcessed(InputWorkspace=sub_result+'_red',Filename=sred_path)
1191                if Verbose:
1192                    logger.notice('Output file created : '+sred_path)
1193            res_plot = sub_result+'_rqw'
1194    if (PlotResult != 'None'):
1195        plotCorrResult(res_plot,PlotResult)
1196    if ( container != '' ):
1197        sws = mtd[sample]
1198        cws = mtd[container]
1199        names = 'Sample,Can,Calc'
1200        for i in range(0, s_hist): # Loop through each spectra in the inputWS
1201            dataX = np.array(sws.readX(i))
1202            dataY = np.array(sws.readY(i))
1203            dataE = np.array(sws.readE(i))
1204            dataX = np.append(dataX,np.array(cws.readX(i)))
1205            dataY = np.append(dataY,np.array(cws.readY(i)))
1206            dataE = np.append(dataE,np.array(cws.readE(i)))
1207            dataX = np.append(dataX,np.array(rws.readX(i)))
1208            dataY = np.append(dataY,np.array(rws.readY(i)))
1209            dataE = np.append(dataE,np.array(rws.readE(i)))
1210            fout = outNm + str(i)
1211            CreateWorkspace(OutputWorkspace=fout, DataX=dataX, DataY=dataY, DataE=dataE,
1212                Nspec=3, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues=names)
1213            if i == 0:
1214                group = fout
1215            else:
1216                group += ',' + fout
1217        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=outNm[:-1])
1218        if PlotContrib:
1219            plotCorrContrib(outNm+'0',[0,1,2])
1220        if Save:
1221            res_path = os.path.join(workdir,outNm[:-1]+'.nxs')
1222            SaveNexusProcessed(InputWorkspace=outNm[:-1],Filename=res_path)
1223            if Verbose:
1224                logger.notice('Output file created : '+res_path)
1225    EndTime('ApplyCorrections')
1226
1227def plotCorrResult(inWS,PlotResult):
1228    nHist = mtd[inWS].getNumberHistograms()
1229    if (PlotResult == 'Spectrum' or PlotResult == 'Both'):
1230        if nHist >= 10:                       #only plot up to 10 hists
1231            nHist = 10
1232        plot_list = []
1233        for i in range(0, nHist):
1234            plot_list.append(i)
1235        res_plot=mp.plotSpectrum(inWS,plot_list)
1236    if (PlotResult == 'Contour' or PlotResult == 'Both'):
1237        if nHist >= 5:                        #needs at least 5 hists for a contour
1238            mp.importMatrixWorkspace(inWS).plotGraph2D()
1239
1240def plotCorrContrib(plot_list,n):
1241        con_plot=mp.plotSpectrum(plot_list,n)