Ticket #8389: IndirectBayes.py

File IndirectBayes.py, 27.6 KB (added by Samuel Jackson, 7 years ago)

Updated version of the script

Line 
1# Bayes routines
2# Fortran programs use fixed length arrays whereas Python has variable lenght lists
3# Input : the Python list is padded to Fortrans length using procedure PadArray
4# Output : the Fortran numpy array is sliced to Python length using dataY = yout[:ny]
5#
6from IndirectImport import *
7if is_supported_f2py_platform():
8    Er      = import_f2py("erange")
9    QLr     = import_f2py("QLres")
10    QLd     = import_f2py("QLdata")
11    Qse     = import_f2py("QLse")
12    Que     = import_f2py("Quest")
13    resnorm = import_f2py("ResNorm")
14else:
15    unsupported_message()
16
17from mantid.simpleapi import *
18from mantid import config, logger, mtd
19from IndirectCommon import *
20import sys, platform, math, os.path, numpy as np
21mp = import_mantidplot()
22
23def CalcErange(inWS,ns,er,nbin):
24        rscl = 1.0
25        array_len = 4096                           # length of array in Fortran
26        N,X,Y,E = GetXYE(inWS,ns,array_len)        # get data
27        nout,bnorm,Xdat=Er.erange(N,X,Y,E,er,nbin,rscl)    # calculate energy range
28        if nout[0] == 0:
29                error = 'Erange - calculated points is Zero'                   
30                logger.notice('ERROR *** ' + error)
31                sys.exit(error)
32        if nout[1] == 0:
33                error = 'Erange - calculated Imin is Zero'                     
34                logger.notice('ERROR *** ' + error)
35                sys.exit(error)
36        if nout[2] == 0:
37                error = 'Erange - calculated Imax is Zero'                     
38                logger.notice('ERROR *** ' + error)
39                sys.exit(error)
40        return nout,bnorm,Xdat,X,Y,E
41
42def GetXYE(inWS,n,array_len):
43        Xin = mtd[inWS].readX(n)
44        N = len(Xin)-1                                                  # get no. points from length of x array
45        Yin = mtd[inWS].readY(n)
46        Ein = mtd[inWS].readE(n)
47        X=PadArray(Xin,array_len)
48        Y=PadArray(Yin,array_len)
49        E=PadArray(Ein,array_len)
50        return N,X,Y,E
51
52def GetResNorm(resnormWS,ngrp):
53        if ngrp == 0:                                # read values from WS
54                dtnorm = mtd[resnormWS+'_Intensity'].readY(0)
55                xscale = mtd[resnormWS+'_Stretch'].readY(0)
56        else:                                        # constant values
57                dtnorm = []
58                xscale = []
59                for m in range(0,ngrp):
60                        dtnorm.append(1.0)
61                        xscale.append(1.0)
62        dtn=PadArray(dtnorm,51)                      # pad for Fortran call
63        xsc=PadArray(xscale,51)
64        return dtn,xsc
65       
66def ReadNormFile(readRes,resnormWS,nsam,Verbose):            # get norm & scale values
67        if readRes:                   # use ResNorm file option=o_res
68                Xin = mtd[resnormWS+'_Intensity'].readX(0)
69                nrm = len(Xin)                                          # no. points from length of x array
70                if nrm == 0:                           
71                        error = 'ResNorm file has no Intensity points'                 
72                        logger.notice('ERROR *** ' + error)
73                        sys.exit(error)
74                Xin = mtd[resnormWS+'_Stretch'].readX(0)                                        # no. points from length of x array
75                if len(Xin) == 0:                               
76                        error = 'ResNorm file has no xscale points'
77                        logger.notice('ERROR *** ' + error)
78                        sys.exit(error)
79                if nrm != nsam:                         # check that no. groups are the same
80                        error = 'ResNorm groups (' +str(nrm) + ') not = Sample (' +str(nsam) +')'                       
81                        logger.notice('ERROR *** ' + error)
82                        sys.exit(error)
83                else:
84                        dtn,xsc = GetResNorm(resnormWS,0)
85        else:
86                # do not use ResNorm file
87                dtn,xsc = GetResNorm(resnormWS,nsam)
88        return dtn,xsc
89
90#Reads in a width ASCII file
91def ReadWidthFile(readWidth,widthFile,numSampleGroups,Verbose):
92        widthY = []
93        widthE = []
94
95        if readWidth:
96
97                if Verbose:     
98                        logger.notice('Width file is ' + widthFile)
99
100                # read ascii based width file
101                try:
102                        wfPath = FileFinder.getFullPath(widthFile)
103                        handle = open(wfPath, 'r')
104                        asc = []
105
106                        for line in handle:
107                                line = line.rstrip()
108                                asc.append(line)
109                        handle.close()
110
111                except Exception, e:
112                        error = 'Failed to read width file'     
113                        logger.notice('ERROR *** ' + error)
114                        sys.exit(error)
115
116                numLines = len(asc)
117               
118                if numLines == 0:
119                        error = 'No groups in width file'       
120                        logger.notice('ERROR *** ' + error)
121                        sys.exit(error)
122               
123                if numLines != numSampleGroups:                         # check that no. groups are the same
124                        error = 'Width groups (' +str(numLines) + ') not = Sample (' +str(numSampleGroups) +')' 
125                        logger.notice('ERROR *** ' + error)
126                        sys.exit(error)
127
128        else: 
129                # no file: just use constant values
130                widthY = np.zeros(numSampleGroups)
131                widthE = np.zeros(numSampleGroups)
132
133        # pad for Fortran call
134        widthY = PadArray(widthY,51)
135        widthE = PadArray(widthE,51)
136
137        return widthY, widthE
138
139# QLines programs
140def QLRun(program,samWS,resWS,resnormWS,erange,nbins,Fit,wfile,Loop,Verbose,Plot,Save):
141        StartTime(program)
142
143        #expand fit options
144        elastic, background, width, resnorm = Fit
145       
146        #convert true/false to 1/0 for fortran
147        o_el = 1 if elastic else 0
148        o_w1 = 1 if width else 0
149        o_res = 1 if resnorm else 0
150
151        #fortran code uses background choices defined using the following numbers
152        if background == 'Sloping':
153                o_bgd = 2
154        elif background == 'Flat':
155                o_bgd = 1
156        elif background == 'Zero':
157                o_bgd = 0
158
159        fitOp = [o_el, o_bgd, o_w1, o_res]
160
161        workdir = getDefaultWorkingDirectory()
162
163        facility = config['default.facility']
164        array_len = 4096                                                   # length of array in Fortran
165        CheckXrange(erange,'Energy')
166
167        nbin,nrbin = nbins[0], nbins[1]
168
169        if Verbose:
170                logger.notice('Sample is ' + samWS)
171                logger.notice('Resolution is ' + resWS)
172
173        if facility == 'ISIS':
174                CheckAnalysers(samWS,resWS,Verbose)
175                efix = getEfixed(samWS)
176                theta,Q = GetThetaQ(samWS)
177
178        nsam,ntc = CheckHistZero(samWS)
179
180        totalNoSam = nsam
181       
182        #check if we're performing a sequential fit
183        if Loop != True:
184                nsam = 1
185
186        nres,ntr = CheckHistZero(resWS)
187       
188        if program == 'QL':
189                if nres == 1:
190                        prog = 'QLr'                                            # res file
191                else:
192                        prog = 'QLd'                                            # data file
193                        CheckHistSame(samWS,'Sample',resWS,'Resolution')
194        elif program == 'QSe':
195                if nres == 1:
196                        prog = 'QSe'                                            # res file
197                else:
198                        error = 'Stretched Exp ONLY works with RES file'
199                        sys.exit(error)
200
201        if Verbose:
202                logger.notice('Version is ' +prog)
203                logger.notice(' Number of spectra = '+str(nsam))
204                logger.notice(' Erange : '+str(erange[0])+' to '+str(erange[1]))
205
206        Wy,We = ReadWidthFile(width,wfile,totalNoSam,Verbose)
207        dtn,xsc = ReadNormFile(resnorm,resnormWS,totalNoSam,Verbose)
208
209        fname = samWS[:-4] + '_'+ prog
210        probWS = fname + '_Prob'
211        fitWS = fname + '_Fit'
212        datWS = fname + '_Data'
213        wrks=os.path.join(workdir, samWS[:-4])
214        if Verbose:
215                logger.notice(' lptfile : '+wrks+'_'+prog+'.lpt')
216        lwrk=len(wrks)
217        wrks.ljust(140,' ')
218        wrkr=resWS
219        wrkr.ljust(140,' ')
220        wrk = [wrks, wrkr]
221
222        # initialise probability list
223        if program == 'QL':
224                prob0 = []
225                prob1 = []
226                prob2 = []
227        xQ = np.array([Q[0]])
228        for m in range(1,nsam):
229                xQ = np.append(xQ,Q[m])
230        xProb = xQ
231        xProb = np.append(xProb,xQ)
232        xProb = np.append(xProb,xQ)
233        eProb = np.zeros(3*nsam)
234
235        group = ''
236        for m in range(0,nsam):
237                if Verbose:
238                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
239                nsp = m+1
240                nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
241                Ndat = nout[0]
242                Imin = nout[1]
243                Imax = nout[2]
244                if prog == 'QLd':
245                        mm = m
246                else:
247                        mm = 0
248                Nb,Xb,Yb,Eb = GetXYE(resWS,mm,array_len)         # get resolution data
249                numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin]
250                rscl = 1.0
251                reals = [efix, theta[m], rscl, bnorm]
252                if prog == 'QLr':
253                        nd,xout,yout,eout,yfit,yprob=QLr.qlres(numb,Xv,Yv,Ev,reals,fitOp,
254                                                                                                   Xdat,Xb,Yb,Wy,We,dtn,xsc,
255                                                                                                   wrks,wrkr,lwrk)
256                        message = ' Log(prob) : '+str(yprob[0])+' '+str(yprob[1])+' '+str(yprob[2])+' '+str(yprob[3])
257                        if Verbose:
258                                logger.notice(message)
259                if prog == 'QLd':
260                        nd,xout,yout,eout,yfit,yprob=QLd.qldata(numb,Xv,Yv,Ev,reals,fitOp,
261                                                                                                        Xdat,Xb,Yb,Eb,Wy,We,
262                                                                                                        wrks,wrkr,lwrk)
263                        message = ' Log(prob) : '+str(yprob[0])+' '+str(yprob[1])+' '+str(yprob[2])+' '+str(yprob[3])
264                        if Verbose:
265                                logger.notice(message)
266                if prog == 'QSe':
267                        nd,xout,yout,eout,yfit,yprob=Qse.qlstexp(numb,Xv,Yv,Ev,reals,fitOp,
268                                                                                                        Xdat,Xb,Yb,Wy,We,dtn,xsc,
269                                                                                                        wrks,wrkr,lwrk)
270                dataX = xout[:nd]
271                dataX = np.append(dataX,2*xout[nd-1]-xout[nd-2])
272                yfit_list = np.split(yfit[:4*nd],4)
273                dataF0 = yfit_list[0]
274                dataF1 = yfit_list[1]
275                if program == 'QL':
276                        dataF2 = yfit_list[2]
277                        dataF3 = yfit_list[3]
278                dataG = np.zeros(nd)
279                datX = dataX
280                datY = yout[:nd]
281                datE = eout[:nd]
282                datX = np.append(datX,dataX)
283                datY = np.append(datY,dataF1[:nd])
284                datE = np.append(datE,dataG)
285                res1 = dataF1[:nd] - yout[:nd]
286                datX = np.append(datX,dataX)
287                datY = np.append(datY,res1)
288                datE = np.append(datE,dataG)
289                nsp = 3
290                names = 'data,fit.1,diff.1'
291                res_plot = [0, 1, 2]
292                if program == 'QL':
293                        datX = np.append(datX,dataX)
294                        datY = np.append(datY,dataF2[:nd])
295                        datE = np.append(datE,dataG)
296                        res2 = dataF2[:nd] - yout[:nd]
297                        datX = np.append(datX,dataX)
298                        datY = np.append(datY,res2)
299                        datE = np.append(datE,dataG)
300                        nsp += 2
301                        names += ',fit.2,diff.2'
302                        res_plot.append(4)
303                        prob0.append(yprob[0])
304                        prob1.append(yprob[1])
305                        prob2.append(yprob[2])
306
307                # create result workspace
308                fitWS = fname+'_Result'
309                fout = fitWS +'_'+ str(m)
310
311                CreateWorkspace(OutputWorkspace=fout, DataX=datX, DataY=datY, DataE=datE,
312                        Nspec=nsp, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues=names)
313               
314                # append workspace to list of results
315                group += fout + ','
316
317       
318        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS)
319
320        if program == 'QL':
321                yPr0 = np.array([prob0[0]])
322                yPr1 = np.array([prob1[0]])
323                yPr2 = np.array([prob2[0]])
324                for m in range(1,nsam):
325                        yPr0 = np.append(yPr0,prob0[m])
326                        yPr1 = np.append(yPr1,prob1[m])
327                        yPr2 = np.append(yPr2,prob2[m])
328                yProb = yPr0
329                yProb = np.append(yProb,yPr1)
330                yProb = np.append(yProb,yPr2)
331                CreateWorkspace(OutputWorkspace=probWS, DataX=xProb, DataY=yProb, DataE=eProb,
332                        Nspec=3, UnitX='MomentumTransfer')
333                outWS = C2Fw(samWS[:-4],fname)
334                if (Plot != 'None'):
335                        QLPlotQL(fname,Plot,res_plot,Loop)
336        if program == 'QSe':
337                outWS = C2Se(fname)
338                if (Plot != 'None'):
339                        QLPlotQSe(fname,Plot,res_plot,Loop)
340
341        #Add some sample logs to the output workspace
342        AddSampleLog(Workspace=outWS, LogName="Fit Program", LogType="String", LogText=prog)
343        AddSampleLog(Workspace=outWS, LogName="Energy min", LogType="Number", LogText=str(erange[0]))
344        AddSampleLog(Workspace=outWS, LogName="Energy max", LogType="Number", LogText=str(erange[1]))
345        AddSampleLog(Workspace=outWS, LogName="Elastic", LogType="String", LogText=str(elastic))
346
347        AddSampleLog(Workspace=outWS, LogName="ResNorm", LogType="String", LogText=str(resnorm))
348        if resnorm:
349                AddSampleLog(Workspace=outWS, LogName="ResNorm file", LogType="String", LogText=resnormWS)
350       
351        AddSampleLog(Workspace=outWS, LogName="Width", LogType="String", LogText=str(width))
352       
353        if width:
354                AddSampleLog(Workspace=outWS, LogName="Width file", LogType="String", LogText=wfile)
355
356        if Save:
357                fit_path = os.path.join(workdir,fitWS+'.nxs')
358                SaveNexusProcessed(InputWorkspace=fitWS, Filename=fit_path)
359                out_path = os.path.join(workdir, outWS+'.nxs')                                  # path name for nxs file
360                SaveNexusProcessed(InputWorkspace=outWS, Filename=out_path)
361                if Verbose:
362                        logger.notice('Output fit file created : ' + fit_path)
363                        logger.notice('Output paramter file created : ' + out_path)
364        EndTime(program)
365
366
367def LorBlock(a,first,nl):                                 #read Ascii block of Integers
368        line1 = a[first]
369        first += 1
370        val = ExtractFloat(a[first])               #Q,AMAX,HWHM,BSCL,GSCL
371        Q = val[0]
372        AMAX = val[1]
373        HWHM = val[2]
374        BSCL = val[3]
375        GSCL = val[4]
376        first += 1
377        val = ExtractFloat(a[first])               #A0,A1,A2,A4
378        int0 = [AMAX*val[0]]
379        bgd1 = BSCL*val[1]
380        bgd2 = BSCL*val[2]
381        zero = GSCL*val[3]
382        first += 1
383        val = ExtractFloat(a[first])                #AI,FWHM first peak
384        fw = [2.*HWHM*val[1]]
385        int = [AMAX*val[0]]
386        if nl >= 2:
387                first += 1
388                val = ExtractFloat(a[first])            #AI,FWHM second peak
389                fw.append(2.*HWHM*val[1])
390                int.append(AMAX*val[0])
391        if nl == 3:
392                first += 1
393                val = ExtractFloat(a[first])            #AI,FWHM third peak
394                fw.append(2.*HWHM*val[1])
395                int.append(AMAX*val[0])
396        first += 1
397        val = ExtractFloat(a[first])                 #SIG0
398        int0.append(val[0])
399        first += 1
400        val = ExtractFloat(a[first])                  #SIGIK
401        int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
402        first += 1
403        val = ExtractFloat(a[first])                  #SIGFK
404        fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
405        if nl >= 2:                                      # second peak
406                first += 1
407                val = ExtractFloat(a[first])                  #SIGIK
408                int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
409                first += 1
410                val = ExtractFloat(a[first])                  #SIGFK
411                fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
412        if nl == 3:                                       # third peak
413                first += 1
414                val = ExtractFloat(a[first])                  #SIGIK
415                int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
416                first += 1
417                val = ExtractFloat(a[first])                  #SIGFK
418                fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
419        first += 1
420        return first,Q,int0,fw,int                                      #values as list
421
422def C2Fw(prog,sname):
423        workdir = config['defaultsave.directory']
424        outWS = sname+'_Workspace'
425        Vaxis = []
426
427        dataX = np.array([])
428        dataY = np.array([])
429        dataE = np.array([])
430
431        nhist = 0
432        for nl in range(1,4):
433                file = sname + '.ql' +str(nl)
434                handle = open(os.path.join(workdir, file), 'r')
435                asc = []
436                for line in handle:
437                        line = line.rstrip()
438                        asc.append(line)
439                handle.close()
440                lasc = len(asc)
441                var = asc[3].split()                                                    #split line on spaces
442                nspec = var[0]
443                ndat = var[1]
444                var = ExtractInt(asc[6])
445                first = 7
446                Xout = []
447
448                ns = int(nspec)
449
450                YData = [[] for i in range(6)]
451                EData = [[] for i in range(6)]
452
453                for m in range(0,ns):
454                        first,Q,i0,fw,it = LorBlock(asc,first,nl)
455                        Xout.append(Q)
456
457                        for i in range(0,nl):
458                                #collect amplitude and width data
459                                YData[i*2].append(fw[i])
460                                YData[i*2+1].append(it[i])
461                                EData[i*2].append(fw[nl+i])
462                                EData[i*2+1].append(it[nl+i])
463
464                nhist += nl*2
465               
466                for i in range(0,nl):
467                        #append amplitude
468                        dataX = np.append(dataX, np.array(Xout))
469                        dataY = np.append(dataY, np.array(YData[i*2+1]))
470                        dataE = np.append(dataE, np.array(EData[i*2+1]))
471                        Vaxis.append('ampl.'+str(nl)+'.'+str(i+1))
472
473                        #append width
474                        dataX = np.append(dataX, np.array(Xout))
475                        dataY = np.append(dataY, np.array(YData[i*2]))
476                        dataE = np.append(dataE, np.array(EData[i*2]))
477                        Vaxis.append('width.'+str(nl)+'.'+str(i+1))
478
479        CreateWorkspace(OutputWorkspace=outWS, DataX=dataX, DataY=dataY, DataE=dataE, Nspec=nhist,
480                UnitX='MomentumTransfer', VerticalAxisUnit='Text', VerticalAxisValues=Vaxis, YUnitLabel='')
481
482        return outWS
483
484def SeBlock(a,first):                                 #read Ascii block of Integers
485        line1 = a[first]
486        first += 1
487        val = ExtractFloat(a[first])               #Q,AMAX,HWHM
488        Q = val[0]
489        AMAX = val[1]
490        HWHM = val[2]
491        first += 1
492        val = ExtractFloat(a[first])               #A0
493        int0 = [AMAX*val[0]]
494        first += 1
495        val = ExtractFloat(a[first])                #AI,FWHM first peak
496        fw = [2.*HWHM*val[1]]
497        int = [AMAX*val[0]]
498        first += 1
499        val = ExtractFloat(a[first])                 #SIG0
500        int0.append(val[0])
501        first += 1
502        val = ExtractFloat(a[first])                  #SIG3K
503        int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
504        first += 1
505        val = ExtractFloat(a[first])                  #SIG1K
506        fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
507        first += 1
508        be = ExtractFloat(a[first])                  #EXPBET
509        first += 1
510        val = ExtractFloat(a[first])                  #SIG2K
511        be.append(math.sqrt(math.fabs(val[0])+1.0e-20))
512        first += 1
513        return first,Q,int0,fw,int,be                                      #values as list
514       
515def C2Se(sname):
516        workdir = config['defaultsave.directory']
517        prog = 'QSe'
518        outWS = sname+'_Workspace'
519        handle = open(os.path.join(workdir, sname+'.qse'), 'r')
520        asc = []
521        for line in handle:
522                line = line.rstrip()
523                asc.append(line)
524        handle.close()
525        lasc = len(asc)
526        var = asc[3].split()                                                    #split line on spaces
527        nspec = var[0]
528        ndat = var[1]
529        var = ExtractInt(asc[6])
530        first = 7
531        Xout = []
532        Yf = []
533        Ef = []
534        Yi = []
535        Ei = []
536        Yb = []
537        Eb = []
538        ns = int(nspec)
539
540        dataX = np.array([])
541        dataY = np.array([])
542        dataE = np.array([])
543
544        for m in range(0,ns):
545                first,Q,int0,fw,it,be = SeBlock(asc,first)
546                Xout.append(Q)
547                Yf.append(fw[0])
548                Ef.append(fw[1])
549                Yi.append(it[0])
550                Ei.append(it[1])
551                Yb.append(be[0])
552                Eb.append(be[1])
553        Vaxis = []
554
555        dataX = np.append(dataX,np.array(Xout))
556        dataY = np.append(dataY,np.array(Yi))
557        dataE = np.append(dataE,np.array(Ei))
558        nhist = 1
559        Vaxis.append('ampl')
560
561        dataX = np.append(dataX, np.array(Xout))
562        dataY = np.append(dataY, np.array(Yf))
563        dataE = np.append(dataE, np.array(Ef))
564        nhist += 1
565        Vaxis.append('width')
566
567        dataX = np.append(dataX,np.array(Xout))
568        dataY = np.append(dataY,np.array(Yb))
569        dataE = np.append(dataE,np.array(Eb))
570        nhist += 1
571        Vaxis.append('beta')
572
573        logger.notice('Vaxis=' + str(Vaxis))
574        CreateWorkspace(OutputWorkspace=outWS, DataX=dataX, DataY=dataY, DataE=dataE, Nspec=nhist,
575                UnitX='MomentumTransfer', VerticalAxisUnit='Text', VerticalAxisValues=Vaxis, YUnitLabel='')
576        return outWS
577
578def QLPlotQL(inputWS,Plot,res_plot,Loop):
579        if Loop:
580                if (Plot == 'Prob' or Plot == 'All'):
581                        pWS = inputWS+'_Prob'
582                        p_plot=mp.plotSpectrum(pWS,[1,2],False)
583
584                if (Plot == 'FwHm' or Plot == 'All'):
585                        ilist = [1,3,5]
586                        i_plot=mp.plotSpectrum(inputWS+'_Workspace',ilist,True)
587                        i_layer = i_plot.activeLayer()
588                        i_layer.setAxisTitle(mp.Layer.Left,'Amplitude')
589                if (Plot == 'Intensity' or Plot == 'All'):
590                        wlist = [0,2,4]
591                        w_plot=mp.plotSpectrum(inputWS+'_Workspace',wlist,True)
592                        w_layer = w_plot.activeLayer()
593                        w_layer.setAxisTitle(mp.Layer.Left,'Full width half maximum (meV)')
594        if (Plot == 'Fit' or Plot == 'All'):
595                fWS = inputWS+'_Result_0'
596                f_plot=mp.plotSpectrum(fWS,res_plot,False)
597
598def QLPlotQSe(inputWS,Plot,res_plot,Loop):
599        if Loop:
600                if (Plot == 'FwHm' or Plot == 'All'):
601                        i_plot=mp.plotSpectrum(inputWS+'_Workspace',1,True)
602                        i_layer = i_plot.activeLayer()
603                        i_layer.setAxisTitle(mp.Layer.Left,'Amplitude')
604                if (Plot == 'Intensity' or Plot == 'All'):
605                        w_plot=mp.plotSpectrum(inputWS+'_Workspace',0,True)
606                        w_layer = w_plot.activeLayer()
607                        w_layer.setAxisTitle(mp.Layer.Left,'Full width half maximum (meV)')
608                if (Plot == 'Beta' or Plot == 'All'):
609                        b_plot=mp.plotSpectrum(inputWS+'_Workspace',2,True)
610                        b_layer = b_plot.activeLayer()
611                        b_layer.setAxisTitle(mp.Layer.Left,'Beta')
612        if (Plot == 'Fit' or Plot == 'All'):
613                fWS = inputWS+'_Result_0'
614                f_plot=mp.plotSpectrum(fWS,res_plot,False)
615
616# Quest programs
617def CheckBetSig(nbs):
618        Nsig = int(nbs[1])
619        if Nsig == 0:
620                error = 'Number of sigma points is Zero'                       
621                logger.notice('ERROR *** ' + error)
622                sys.exit(error)
623        if Nsig > 200:
624                error = 'Max number of sigma points is 200'                     
625                logger.notice('ERROR *** ' + error)
626                sys.exit(error)
627        Nbet = int(nbs[0])
628        if Nbet == 0:
629                error = 'Number of beta points is Zero'                 
630                logger.notice('ERROR *** ' + error)
631                sys.exit(error)
632        if Nbet > 200:
633                error = 'Max number of beta points is 200'                     
634                logger.notice('ERROR *** ' + error)
635                sys.exit(error)
636        return Nbet,Nsig
637
638def QuestRun(samWS,resWS,nbs,erange,nbins,Fit,Loop,Verbose,Plot,Save):
639        StartTime('Quest')
640        #expand fit options
641        elastic, background, width, resnorm = Fit
642       
643        #convert true/false to 1/0 for fortran
644        o_el = 1 if elastic else 0
645        o_w1 = 1 if width else 0
646        o_res = 1 if resnorm else 0
647
648        #fortran code uses background choices defined using the following numbers
649        if background == 'Sloping':
650                o_bgd = 2
651        elif background == 'Flat':
652                o_bgd = 1
653        elif background == 'Zero':
654                o_bgd = 0
655
656        fitOp = [o_el, o_bgd, o_w1, o_res]
657
658        workdir = getDefaultWorkingDirectory()
659
660        array_len = 4096                           # length of array in Fortran
661        CheckXrange(erange,'Energy')
662        nbin,nrbin = nbins[0],nbins[1]
663        if Verbose:
664                logger.notice('Sample is ' + samWS)
665                logger.notice('Resolution is ' + resWS)
666        CheckAnalysers(samWS,resWS,Verbose)
667        nsam,ntc = CheckHistZero(samWS)
668       
669        if Loop != True:
670                nsam = 1
671
672        efix = getEfixed(samWS)
673        theta,Q = GetThetaQ(samWS)
674        nres,ntr = CheckHistZero(resWS)
675        if nres == 1:
676                prog = 'Qst'                        # res file
677        else:
678                error = 'Stretched Exp ONLY works with RES file'                       
679                logger.notice('ERROR *** ' + error)
680                sys.exit(error)
681        if Verbose:
682                logger.notice(' Number of spectra = '+str(nsam))
683                logger.notice(' Erange : '+str(erange[0])+' to '+str(erange[1]))
684
685        fname = samWS[:-4] + '_'+ prog
686        wrks=os.path.join(workdir, samWS[:-4])
687        if Verbose:
688                logger.notice(' lptfile : ' + wrks +'_Qst.lpt')
689        lwrk=len(wrks)
690        wrks.ljust(140,' ')
691        wrkr=resWS
692        wrkr.ljust(140,' ')
693        wrk = [wrks, wrkr]
694        Nbet,Nsig = nbs[0], nbs[1]
695        eBet0 = np.zeros(Nbet)                  # set errors to zero
696        eSig0 = np.zeros(Nsig)                  # set errors to zero
697        rscl = 1.0
698        Qaxis = ''
699        for m in range(0,nsam):
700                if Verbose:
701                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
702                nsp = m+1
703                nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
704                Ndat = nout[0]
705                Imin = nout[1]
706                Imax = nout[2]
707                Nb,Xb,Yb,Eb = GetXYE(resWS,0,array_len)
708                numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin, Nbet, Nsig]
709                reals = [efix, theta[m], rscl, bnorm]
710                xsout,ysout,xbout,ybout,zpout=Que.quest(numb,Xv,Yv,Ev,reals,fitOp,
711                                                                                        Xdat,Xb,Yb,wrks,wrkr,lwrk)
712                dataXs = xsout[:Nsig]               # reduce from fixed Fortran array
713                dataYs = ysout[:Nsig]
714                dataXb = xbout[:Nbet]
715                dataYb = ybout[:Nbet]
716                zpWS = fname + '_Zp' +str(m)
717                if (m > 0):
718                        Qaxis += ','
719                Qaxis += str(Q[m])
720                for n in range(0,Nsig):
721                        yfit_list = np.split(zpout[:Nsig*Nbet],Nsig)
722                        dataYzp = yfit_list[n]
723                        if n == 0:
724                                dataXz = xbout[:Nbet]
725                                dataYz = dataYzp[:Nbet]
726                                dataEz = eBet0
727                        else:
728                                dataXz = np.append(dataXz,xbout[:Nbet])
729                                dataYz = np.append(dataYz,dataYzp[:Nbet])
730                                dataEz = np.append(dataEz,eBet0)
731                CreateWorkspace(OutputWorkspace=zpWS, DataX=dataXz, DataY=dataYz, DataE=dataEz,
732                        Nspec=Nsig, UnitX='MomentumTransfer', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=dataXs)
733                unitx = mtd[zpWS].getAxis(0).setUnit("Label")
734                unitx.setLabel('beta' , '')
735                unity = mtd[zpWS].getAxis(1).setUnit("Label")
736                unity.setLabel('sigma' , '')
737                if m == 0:
738                        xSig = dataXs
739                        ySig = dataYs
740                        eSig = eSig0           
741                        xBet = dataXb
742                        yBet = dataYb
743                        eBet = eBet0           
744                        groupZ = zpWS
745                else:
746                        xSig = np.append(xSig,dataXs)
747                        ySig = np.append(ySig,dataYs)
748                        eSig = np.append(eSig,eSig0)
749                        xBet = np.append(xBet,dataXb)
750                        yBet = np.append(yBet,dataYb)
751                        eBet = np.append(eBet,eBet0)
752                        groupZ = groupZ +','+ zpWS
753        CreateWorkspace(OutputWorkspace=fname+'_Sigma', DataX=xSig, DataY=ySig, DataE=eSig,
754                Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
755        unitx = mtd[fname+'_Sigma'].getAxis(0).setUnit("Label")
756        unitx.setLabel('sigma' , '')
757        CreateWorkspace(OutputWorkspace=fname+'_Beta', DataX=xBet, DataY=yBet, DataE=eBet,
758                Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
759        unitx = mtd[fname+'_Beta'].getAxis(0).setUnit("Label")
760        unitx.setLabel('beta' , '')
761        group = fname + '_Sigma,'+ fname + '_Beta'
762        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_Fit')     
763        GroupWorkspaces(InputWorkspaces=groupZ,OutputWorkspace=fname+'_Contour')
764
765        if Save:
766                fpath = os.path.join(workdir,fname+'_Fit.nxs')
767                SaveNexusProcessed(InputWorkspace=fname+'_Fit', Filename=fpath)
768                cpath = os.path.join(workdir,fname+'_Contour.nxs')
769                SaveNexusProcessed(InputWorkspace=fname+'_Contour', Filename=cpath)
770                if Verbose:
771                        logger.notice('Output file for Fit : ' + fpath)
772                        logger.notice('Output file for Contours : ' + cpath)
773        if (Plot != 'None'):
774                QuestPlot(fname,Plot)
775        EndTime('Quest')
776
777def QuestPlot(inputWS,Plot):
778        if (Plot == 'Sigma' or Plot == 'All'):
779                sig_plot=mp.importMatrixWorkspace(inputWS+'_Sigma').plotGraph2D()
780        if (Plot == 'Beta' or Plot == 'All'):
781                beta_plot=mp.importMatrixWorkspace(inputWS+'_Beta').plotGraph2D()
782
783# ResNorm programs
784def ResNormRun(vname,rname,erange,nbin,Verbose=False,Plot='None',Save=False):
785        StartTime('ResNorm')
786       
787        workdir = getDefaultWorkingDirectory()
788
789        array_len = 4096                                    # length of Fortran array
790        CheckXrange(erange,'Energy')
791        CheckAnalysers(vname,rname,Verbose)
792        nvan,ntc = CheckHistZero(vname)
793        theta,Q = GetThetaQ(vname)
794        efix = getEfixed(vname)
795        nres,ntr = CheckHistZero(rname)
796        nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(vname,0,erange,nbin)
797        Ndat = nout[0]
798        Imin = nout[1]
799        Imax = nout[2]
800        wrks=os.path.join(workdir, vname[:-4])
801        if Verbose:
802                logger.notice(' Number of spectra = '+str(nvan))
803                logger.notice(' lptfile : ' + wrks +'_resnrm.lpt')
804        lwrk=len(wrks)
805        wrks.ljust(140,' ')                              # pad for fioxed Fortran length
806        wrkr=rname
807        wrkr.ljust(140,' ')
808        Nb,Xb,Yb,Eb = GetXYE(rname,0,array_len)
809        rscl = 1.0
810        xPar = np.array([theta[0]])
811        for m in range(1,nvan):
812                xPar = np.append(xPar,theta[m])
813        ePar = np.zeros(nvan)
814        fname = vname[:-4]
815        for m in range(0,nvan):
816                if Verbose:
817                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
818                ntc,Xv,Yv,Ev = GetXYE(vname,m,array_len)
819                nsp = m+1
820                numb = [nvan, nsp, ntc, Ndat, nbin, Imin, Imax, Nb]
821                reals = [efix, theta[0], rscl, bnorm]
822                nd,xout,yout,eout,yfit,pfit=resnorm.resnorm(numb,Xv,Yv,Ev,reals,
823                                                                        Xdat,Xb,Yb,wrks,wrkr,lwrk)
824                if Verbose:
825                        message = ' Fit paras : '+str(pfit[0])+' '+str(pfit[1])
826                        logger.notice(message)
827                dataX = xout[:nd]
828                dataX = np.append(dataX,2*xout[nd-1]-xout[nd-2])
829                if m == 0:
830                        yPar1 = np.array([pfit[0]])
831                        yPar2 = np.array([pfit[1]])
832                        CreateWorkspace(OutputWorkspace='Data', DataX=dataX, DataY=yout[:nd], DataE=eout[:nd],
833                                NSpec=1, UnitX='DeltaE')
834                        CreateWorkspace(OutputWorkspace='Fit', DataX=dataX, DataY=yfit[:nd], DataE=np.zeros(nd),
835                                NSpec=1, UnitX='DeltaE')
836                else:
837                        yPar1 = np.append(yPar1,pfit[0])
838                        yPar2 = np.append(yPar2,pfit[1])
839                        CreateWorkspace(OutputWorkspace='__datmp', DataX=dataX, DataY=yout[:nd], DataE=eout[:nd],
840                                NSpec=1, UnitX='DeltaE')
841                        ConjoinWorkspaces(InputWorkspace1='Data', InputWorkspace2='__datmp', CheckOverlapping=False)                           
842                        CreateWorkspace(OutputWorkspace='__f1tmp', DataX=dataX, DataY=yfit[:nd], DataE=np.zeros(nd),
843                                NSpec=1, UnitX='DeltaE')
844                        ConjoinWorkspaces(InputWorkspace1='Fit', InputWorkspace2='__f1tmp', CheckOverlapping=False)                             
845        CreateWorkspace(OutputWorkspace=fname+'_ResNorm_Intensity', DataX=xPar, DataY=yPar1, DataE=xPar,
846                NSpec=1, UnitX='MomentumTransfer')
847        CreateWorkspace(OutputWorkspace=fname+'_ResNorm_Stretch', DataX=xPar, DataY=yPar2, DataE=xPar,
848                NSpec=1, UnitX='MomentumTransfer')
849        group = fname + '_ResNorm_Intensity,'+ fname + '_ResNorm_Stretch'
850        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_ResNorm')
851        GroupWorkspaces(InputWorkspaces='Data,Fit',OutputWorkspace=fname+'_ResNorm_Fit')
852        if Save:
853                par_path = os.path.join(workdir,fname+'_ResNorm.nxs')
854                SaveNexusProcessed(InputWorkspace=fname+'_ResNorm', Filename=par_path)
855                fit_path = os.path.join(workdir,fname+'_ResNorm_Fit.nxs')
856                SaveNexusProcessed(InputWorkspace=fname+'_ResNorm_Fit', Filename=fit_path)
857                if Verbose:
858                        logger.notice('Parameter file created : ' + par_path)
859                        logger.notice('Fit file created : ' + fit_path)
860        if (Plot != 'None'):
861                ResNormPlot(fname,Plot)
862        EndTime('ResNorm')
863
864def ResNormPlot(inputWS,Plot):
865        if (Plot == 'Intensity' or Plot == 'All'):
866                iWS = inputWS + '_ResNorm_Intensity'
867                i_plot=mp.plotSpectrum(iWS,0,False)
868        if (Plot == 'Stretch' or Plot == 'All'):
869                sWS = inputWS + '_ResNorm_Stretch'
870                s_plot=mp.plotSpectrum(sWS,0,False)
871        if (Plot == 'Fit' or Plot == 'All'):
872                fWS = inputWS + '_ResNorm_Fit'
873                f_plot=mp.plotSpectrum(fWS,0,False)