Ticket #7966: IndirectBayes.py

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

Updated version of IndirectBayes

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(ngrp):
53        if ngrp == 0:                                # read values from WS
54                dtnorm = mtd['ResNorm_1'].readY(0)
55                xscale = mtd['ResNorm_2'].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(o_res,nsam,Verbose):            # get norm & scale values
67        if o_res == 1:                     # use ResNorm file option=o_res
68                Xin = mtd['ResNorm_1'].readX(0)
69                nrm = len(Xin)                                          # no. points from length of x array
70                if nrm == 0:                           
71                        error = 'ResNorm file has no dtnorm points'                     
72                        logger.notice('ERROR *** ' + error)
73                        sys.exit(error)
74                Xin = mtd['ResNorm_2'].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(0)
85        if o_res == 0:                     # do not use ResNorm file option=o_res
86                dtn,xsc = GetResNorm(nsam)
87        return dtn,xsc
88
89def ReadWidthFile(op_w1,wfile,ngrp,Verbose):                       # reads width file ASCII
90        workdir = config['defaultsave.directory']
91        if op_w1 == 1:                               # use width1 data  option=o_w1
92                if Verbose:
93                        w_path = os.path.join(workdir, wfile)                                   # path name for nxs file
94                        logger.notice('Width file is ' + w_path)
95                handle = open(w_path, 'r')
96                asc = []
97                for line in handle:
98                        line = line.rstrip()
99                        asc.append(line)
100                handle.close()
101                lasc = len(asc)
102                if lasc == 0:
103                        error = 'No groups in width file'       
104                        logger.notice('ERROR *** ' + error)
105                        sys.exit(error)
106                if lasc != ngrp:                                # check that no. groups are the same
107                        error = 'Width groups (' +str(lasc) + ') not = Sample (' +str(ngrp) +')'       
108                        logger.notice('ERROR *** ' + error)
109                        sys.exit(error)
110        else:                                           # constant values
111                Wy = []
112                We = []
113                for m in range(0,ngrp):
114                        Wy.append(0.0)
115                        We.append(0.0)
116        Wy=PadArray(Wy,51)                             # pad for Fortran call
117        We=PadArray(We,51)
118        return Wy,We
119
120def CheckBinning(nbins):
121        nbin = nbins[0]
122        if nbin == '0':
123                error = 'Sample binning is Zero'                       
124                logger.notice('ERROR *** ' + error)
125                sys.exit(error)
126        if len(nbins) == 2:
127                nrbin = nbins[1]
128                if nrbin == '0':
129                        error = 'Resolution binning is Zero'                   
130                        logger.notice('ERROR *** ' + error)
131                        sys.exit(error)
132        else:
133                nrbin = 1
134        return nbin,nrbin
135
136# QLines programs
137       
138def QLRun(program,samWS,resWS,rsname,erange,nbins,fitOp,wfile,Loop,Verbose,Plot,Save):
139        StartTime(program)
140        workdir = config['defaultsave.directory']
141        facility = config['default.facility']
142        array_len = 4096                           # length of array in Fortran
143        CheckXrange(erange,'Energy')
144        nbin,nrbin = CheckBinning(nbins)
145        if Verbose:
146                logger.notice('Sample is ' + samWS)
147                logger.notice('Resolution is ' + resWS)
148        if facility == 'ISIS':
149                CheckAnalysers(samWS,resWS,Verbose)
150                efix = getEfixed(samWS)
151                theta,Q = GetThetaQ(samWS)
152        nsam,ntc = CheckHistZero(samWS)
153        if Loop != True:
154                nsam = 1
155        nres,ntr = CheckHistZero(resWS)
156        if program == 'QL':
157                if nres == 1:
158                        prog = 'QLr'                        # res file
159                else:
160                        prog = 'QLd'                        # data file
161                        CheckHistSame(samWS,'Sample',resWS,'Resolution')
162        if program == 'QSe':
163                if nres == 1:
164                        prog = 'QSe'                        # res file
165                else:
166                        error = 'Stretched Exp ONLY works with RES file'                       
167                        logger.notice('ERROR *** ' + error)
168                        sys.exit(error)
169        if Verbose:
170                logger.notice('Version is ' +prog)
171                logger.notice(' Number of spectra = '+str(nsam))
172                logger.notice(' Erange : '+str(erange[0])+' to '+str(erange[1]))
173        Wy,We = ReadWidthFile(fitOp[2],wfile,nsam,Verbose)
174        dtn,xsc = ReadNormFile(fitOp[3],nsam,Verbose)
175        fname = samWS[:-4] + '_'+ prog
176        probWS = fname + '_Prob'
177        fitWS = fname + '_Fit'
178        datWS = fname + '_Data'
179        if program == 'QSe':
180                fwWS = fname + '_FwHm'
181                itWS = fname + '_Inty'
182                beWS = fname + '_Beta'
183        if program == 'QL':
184                fit1WS = fname + '_Fit1'
185                fit2WS = fname + '_Fit2'
186                res1WS = fname + '_Res1'
187                res2WS = fname + '_Res2'
188        wrks=workdir + samWS[:-4]
189        if Verbose:
190                logger.notice(' lptfile : '+wrks+'_'+prog+'.lpt')
191        lwrk=len(wrks)
192        wrks.ljust(140,' ')
193        wrkr=resWS
194        wrkr.ljust(140,' ')
195        wrk = [wrks, wrkr]
196#
197        if program == 'QL':                           # initialise probability list
198                prob0 = []
199                prob1 = []
200                prob2 = []
201        xQ = np.array([Q[0]])
202        for m in range(1,nsam):
203                xQ = np.append(xQ,Q[m])
204        xProb = xQ
205        xProb = np.append(xProb,xQ)
206        xProb = np.append(xProb,xQ)
207        eProb = np.zeros(3*nsam)
208        for m in range(0,nsam):
209                if Verbose:
210                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
211                nsp = m+1
212                nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
213                Ndat = nout[0]
214                Imin = nout[1]
215                Imax = nout[2]
216                if prog == 'QLd':
217                        mm = m
218                else:
219                        mm = 0
220                Nb,Xb,Yb,Eb = GetXYE(resWS,mm,array_len)     # get resolution data
221                numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin]
222                rscl = 1.0
223                reals = [efix, theta[m], rscl, bnorm]
224                if prog == 'QLr':
225                        nd,xout,yout,eout,yfit,yprob=QLr.qlres(numb,Xv,Yv,Ev,reals,fitOp,
226                                                                                                   Xdat,Xb,Yb,Wy,We,dtn,xsc,
227                                                                                                   wrks,wrkr,lwrk)
228                        message = ' Log(prob) : '+str(yprob[0])+' '+str(yprob[1])+' '+str(yprob[2])+' '+str(yprob[3])
229                        if Verbose:
230                                logger.notice(message)
231                if prog == 'QLd':
232                        nd,xout,yout,eout,yfit,yprob=QLd.qldata(numb,Xv,Yv,Ev,reals,fitOp,
233                                                                                                    Xdat,Xb,Yb,Eb,Wy,We,
234                                                                                                        wrks,wrkr,lwrk)
235                        message = ' Log(prob) : '+str(yprob[0])+' '+str(yprob[1])+' '+str(yprob[2])+' '+str(yprob[3])
236                        if Verbose:
237                                logger.notice(message)
238                if prog == 'QSe':
239                        nd,xout,yout,eout,yfit,yprob=Qse.qlstexp(numb,Xv,Yv,Ev,reals,fitOp,
240                                                                                                        Xdat,Xb,Yb,Wy,We,dtn,xsc,
241                                                                                                        wrks,wrkr,lwrk)
242                dataX = xout[:nd]
243                dataX = np.append(dataX,2*xout[nd-1]-xout[nd-2])
244                yfit_list = np.split(yfit[:4*nd],4)
245                dataF0 = yfit_list[0]
246                dataF1 = yfit_list[1]
247                if program == 'QL':
248                        dataF2 = yfit_list[2]
249                        dataF3 = yfit_list[3]
250                dataG = np.zeros(nd)
251                datX = dataX
252                datY = yout[:nd]
253                datE = eout[:nd]
254                datX = np.append(datX,dataX)
255                datY = np.append(datY,dataF1[:nd])
256                datE = np.append(datE,dataG)
257                res1 = dataF1[:nd] - yout[:nd]
258                datX = np.append(datX,dataX)
259                datY = np.append(datY,res1)
260                datE = np.append(datE,dataG)
261                nsp = 3
262                names = 'data,fit.1,diff.1'
263                res_plot = [0, 1, 2]
264                if program == 'QL':
265                        datX = np.append(datX,dataX)
266                        datY = np.append(datY,dataF2[:nd])
267                        datE = np.append(datE,dataG)
268                        res2 = dataF2[:nd] - yout[:nd]
269                        datX = np.append(datX,dataX)
270                        datY = np.append(datY,res2)
271                        datE = np.append(datE,dataG)
272                        nsp += 2
273                        names += ',fit.2,diff.2'
274                        res_plot.append(4)
275                        prob0.append(yprob[0])
276                        prob1.append(yprob[1])
277                        prob2.append(yprob[2])
278                fitWS = fname+'_Result'
279                fout = fitWS +'_'+ str(m)
280                CreateWorkspace(OutputWorkspace=fout, DataX=datX, DataY=datY, DataE=datE,
281                        Nspec=nsp, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues=names)
282                if m == 0:
283                        group = fout
284                else:
285                        group += ',' + fout
286        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS)
287        if program == 'QL':
288                yPr0 = np.array([prob0[0]])
289                yPr1 = np.array([prob1[0]])
290                yPr2 = np.array([prob2[0]])
291                for m in range(1,nsam):
292                        yPr0 = np.append(yPr0,prob0[m])
293                        yPr1 = np.append(yPr1,prob1[m])
294                        yPr2 = np.append(yPr2,prob2[m])
295                yProb = yPr0
296                yProb = np.append(yProb,yPr1)
297                yProb = np.append(yProb,yPr2)
298                CreateWorkspace(OutputWorkspace=probWS, DataX=xProb, DataY=yProb, DataE=eProb,
299                        Nspec=3, UnitX='MomentumTransfer')
300                C2Fw(samWS[:-4],fname)
301                if (Plot != 'None'):
302                        QLPlotQL(fname,Plot,res_plot,Loop)
303        if program == 'QSe':
304                C2Se(fname)
305                if (Plot != 'None'):
306                        QLPlotQSe(fname,Plot,res_plot,Loop)
307        if Save:
308                fit_path = os.path.join(workdir,fitWS+'.nxs')
309                SaveNexusProcessed(InputWorkspace=fitWS, Filename=fit_path)
310                if Verbose:
311                        logger.notice('Output file created : ' + fit_path)
312        EndTime(program)
313
314def LorBlock(a,first,nl):                                 #read Ascii block of Integers
315        line1 = a[first]
316        first += 1
317        val = ExtractFloat(a[first])               #Q,AMAX,HWHM,BSCL,GSCL
318        Q = val[0]
319        AMAX = val[1]
320        HWHM = val[2]
321        BSCL = val[3]
322        GSCL = val[4]
323        first += 1
324        val = ExtractFloat(a[first])               #A0,A1,A2,A4
325        int0 = [AMAX*val[0]]
326        bgd1 = BSCL*val[1]
327        bgd2 = BSCL*val[2]
328        zero = GSCL*val[3]
329        first += 1
330        val = ExtractFloat(a[first])                #AI,FWHM first peak
331        fw = [2.*HWHM*val[1]]
332        int = [AMAX*val[0]]
333        if nl >= 2:
334                first += 1
335                val = ExtractFloat(a[first])            #AI,FWHM second peak
336                fw.append(2.*HWHM*val[1])
337                int.append(AMAX*val[0])
338        if nl == 3:
339                first += 1
340                val = ExtractFloat(a[first])            #AI,FWHM third peak
341                fw.append(2.*HWHM*val[1])
342                int.append(AMAX*val[0])
343        first += 1
344        val = ExtractFloat(a[first])                 #SIG0
345        int0.append(val[0])
346        first += 1
347        val = ExtractFloat(a[first])                  #SIGIK
348        int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
349        first += 1
350        val = ExtractFloat(a[first])                  #SIGFK
351        fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
352        if nl >= 2:                                      # second peak
353                first += 1
354                val = ExtractFloat(a[first])                  #SIGIK
355                int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
356                first += 1
357                val = ExtractFloat(a[first])                  #SIGFK
358                fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
359        if nl == 3:                                       # third peak
360                first += 1
361                val = ExtractFloat(a[first])                  #SIGIK
362                int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
363                first += 1
364                val = ExtractFloat(a[first])                  #SIGFK
365                fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
366        first += 1
367        return first,Q,int0,fw,int                                      #values as list
368       
369def ReadQlFile(prog,sname,nl):
370        workdir = config['defaultsave.directory']
371        fname = sname
372        file = fname + '.ql' +str(nl)
373        handle = open(os.path.join(workdir, file), 'r')
374        asc = []
375        for line in handle:
376                line = line.rstrip()
377                asc.append(line)
378        handle.close()
379        lasc = len(asc)
380        var = asc[3].split()                                                    #split line on spaces
381        nspec = var[0]
382        ndat = var[1]
383        var = ExtractInt(asc[6])
384        first = 7
385        Xout = []
386        Yf1 = []
387        Ef1 = []
388        Yf2 = []
389        Ef2 = []
390        Yf3 = []
391        Ef3 = []
392        Yi1 = []
393        Ei1 = []
394        Yi2 = []
395        Ei2 = []
396        Yi3 = []
397        Ei3 = []
398        ns = int(nspec)
399        for m in range(0,ns):
400                if nl == 1:
401                        first,Q,int0,fw,it = LorBlock(asc,first,1)
402                        Xout.append(Q)
403                        Yf1.append(fw[0])
404                        Ef1.append(fw[1])
405                        Yi1.append(it[0])
406                        Ei1.append(it[1])
407                if nl == 2:
408                        first,Q,int0,fw,it = LorBlock(asc,first,2)
409                        Xout.append(Q)
410                        Yf1.append(fw[0])
411                        Ef1.append(fw[2])
412                        Yf2.append(fw[1])
413                        Ef2.append(fw[3])
414                        Yi1.append(it[0])
415                        Ei1.append(it[2])
416                        Yi2.append(it[1])
417                        Ei2.append(it[3])
418                if nl == 3:
419                        first,Q,int0,fw,it = LorBlock(asc,first,3)
420                        Xout.append(Q)
421                        Yf1.append(fw[0])
422                        Ef1.append(fw[3])
423                        Yf2.append(fw[1])
424                        Ef2.append(fw[4])
425                        Yf3.append(fw[2])
426                        Ef3.append(fw[5])
427                        Yi1.append(it[0])
428                        Ei1.append(it[3])
429                        Yi2.append(it[1])
430                        Ei2.append(it[4])
431                        Yi3.append(it[2])
432                        Ei3.append(it[5])
433        if nl == 1:
434                CreateWorkspace(OutputWorkspace=fname+'_FW11', DataX=Xout, DataY=Yf1, DataE=Ef1,
435                        Nspec=1, UnitX='MomentumTransfer')
436                CreateWorkspace(OutputWorkspace=fname+'_IT11', DataX=Xout, DataY=Yi1, DataE=Ei1,
437                        Nspec=1, UnitX='MomentumTransfer')
438        if nl == 2:
439                CreateWorkspace(OutputWorkspace=fname+'_FW21', DataX=Xout, DataY=Yf1, DataE=Ef1,
440                        Nspec=1, UnitX='MomentumTransfer')
441                CreateWorkspace(OutputWorkspace=fname+'_FW22', DataX=Xout, DataY=Yf2, DataE=Ef2,
442                        Nspec=1, UnitX='MomentumTransfer')
443                CreateWorkspace(OutputWorkspace=fname+'_IT21', DataX=Xout, DataY=Yi1, DataE=Ei1,
444                        Nspec=1, UnitX='MomentumTransfer')
445                CreateWorkspace(OutputWorkspace=fname+'_IT22', DataX=Xout, DataY=Yi2, DataE=Ei2,
446                        Nspec=1, UnitX='MomentumTransfer')
447        if nl == 3:
448                CreateWorkspace(OutputWorkspace=fname+'_FW31', DataX=Xout, DataY=Yf1, DataE=Ef1,
449                        Nspec=1, UnitX='MomentumTransfer')
450                CreateWorkspace(OutputWorkspace=fname+'_FW32', DataX=Xout, DataY=Yf2, DataE=Ef2,
451                        Nspec=1, UnitX='MomentumTransfer')
452                CreateWorkspace(OutputWorkspace=fname+'_FW33', DataX=Xout, DataY=Yf3, DataE=Ef3,       
453                        Nspec=1, UnitX='MomentumTransfer')
454                CreateWorkspace(OutputWorkspace=fname+'_IT31', DataX=Xout, DataY=Yi1, DataE=Ei1,       
455                        Nspec=1, UnitX='MomentumTransfer')
456                CreateWorkspace(OutputWorkspace=fname+'_IT32', DataX=Xout, DataY=Yi2, DataE=Ei2,
457                        Nspec=1, UnitX='MomentumTransfer')
458                CreateWorkspace(OutputWorkspace=fname+'_IT33', DataX=Xout, DataY=Yi3, DataE=Ei3,
459                        Nspec=1, UnitX='MomentumTransfer')
460
461def C2Fw(prog,sname):
462        workdir = config['defaultsave.directory']
463        fname = sname
464        ReadQlFile(prog,sname,1)
465        ReadQlFile(prog,sname,2)
466        ReadQlFile(prog,sname,3)
467        group2 = fname + '_FW21,'+ fname + '_FW22'
468        group3 = fname + '_FW31,'+ fname + '_FW32,'+ fname + '_FW33'
469        group = fname + '_FW11,'+  group2 +','+  group3
470        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_FwHm')
471        group2 = fname + '_IT21,'+ fname + '_IT22'
472        group3 = fname + '_IT31,'+ fname + '_IT32,'+ fname + '_IT33'
473        group = fname + '_IT11,'+ group2 +','+ group3
474        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_Inty')
475        group = fname + '_FwHm,'+ fname + '_Inty'
476        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_Parameters')
477        opath = os.path.join(workdir,fname+'_Parameters.nxs')
478        SaveNexusProcessed(InputWorkspace=fname+'_Parameters', Filename=opath)
479
480def SeBlock(a,first):                                 #read Ascii block of Integers
481        line1 = a[first]
482        first += 1
483        val = ExtractFloat(a[first])               #Q,AMAX,HWHM
484        Q = val[0]
485        AMAX = val[1]
486        HWHM = val[2]
487        first += 1
488        val = ExtractFloat(a[first])               #A0
489        int0 = [AMAX*val[0]]
490        first += 1
491        val = ExtractFloat(a[first])                #AI,FWHM first peak
492        fw = [2.*HWHM*val[1]]
493        int = [AMAX*val[0]]
494        first += 1
495        val = ExtractFloat(a[first])                 #SIG0
496        int0.append(val[0])
497        first += 1
498        val = ExtractFloat(a[first])                  #SIG3K
499        int.append(AMAX*math.sqrt(math.fabs(val[0])+1.0e-20))
500        first += 1
501        val = ExtractFloat(a[first])                  #SIG1K
502        fw.append(2.0*HWHM*math.sqrt(math.fabs(val[0])+1.0e-20))
503        first += 1
504        be = ExtractFloat(a[first])                  #EXPBET
505        first += 1
506        val = ExtractFloat(a[first])                  #SIG2K
507        be.append(math.sqrt(math.fabs(val[0])+1.0e-20))
508        first += 1
509        return first,Q,int0,fw,int,be                                      #values as list
510       
511def C2Se(sname):
512        workdir = config['defaultsave.directory']
513        prog = 'QSe'
514        fname = sname
515        file = fname + '.qse'
516        handle = open(os.path.join(workdir, file), 'r')
517        asc = []
518        for line in handle:
519                line = line.rstrip()
520                asc.append(line)
521        handle.close()
522        lasc = len(asc)
523        var = asc[3].split()                                                    #split line on spaces
524        nspec = var[0]
525        ndat = var[1]
526        var = ExtractInt(asc[6])
527        first = 7
528        Xout = []
529        Yf = []
530        Ef = []
531        Yi = []
532        Ei = []
533        Yb = []
534        Eb = []
535        ns = int(nspec)
536        for m in range(0,ns):
537                first,Q,int0,fw,it,be = SeBlock(asc,first)
538                Xout.append(Q)
539                Yf.append(fw[0])
540                Ef.append(fw[1])
541                Yi.append(it[0])
542                Ei.append(it[1])
543                Yb.append(be[0])
544                Eb.append(be[1])
545        CreateWorkspace(OutputWorkspace=fname+'_FwHm', DataX=Xout, DataY=Yf, DataE=Ef,
546                Nspec=1, UnitX='MomentumTransfer')
547        CreateWorkspace(OutputWorkspace=fname+'_Inty', DataX=Xout, DataY=Yi, DataE=Ei,
548                Nspec=1, UnitX='MomentumTransfer')
549        CreateWorkspace(OutputWorkspace=fname+'_Beta', DataX=Xout, DataY=Yb, DataE=Eb,
550                Nspec=1, UnitX='MomentumTransfer')
551        group = fname + '_FwHm,'+ fname + '_Inty,'+ fname + '_Beta'
552        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_Parameters')
553        opath = os.path.join(workdir,fname+'_Parameters.nxs')
554        SaveNexusProcessed(InputWorkspace=fname+'_Parameters', Filename=opath)
555
556def QLPlotQL(inputWS,Plot,res_plot,Loop):
557        if Loop:
558                if (Plot == 'Prob' or Plot == 'All'):
559                        pWS = inputWS+'_Prob'
560                        p_plot=mp.plotSpectrum(pWS,[1,2],False)
561                if (Plot == 'Intensity' or Plot == 'All'):
562                        iWS = [inputWS+'_IT11', inputWS+'_IT21', inputWS+'_IT22']
563                        i_plot=mp.plotSpectrum(iWS,0,True)
564                if (Plot == 'FwHm' or Plot == 'All'):
565                        wWS = [inputWS+'_FW11', inputWS+'_FW21', inputWS+'_FW22']
566                        w_plot=mp.plotSpectrum(wWS,0,True)
567        if (Plot == 'Fit' or Plot == 'All'):
568                fWS = inputWS+'_Result_0'
569                f_plot=mp.plotSpectrum(fWS,res_plot,False)
570
571def QLPlotQSe(inputWS,Plot,res_plot,Loop):
572        if Loop:
573                if (Plot == 'Intensity' or Plot == 'All'):
574                        i_plot=mp.plotSpectrum(inputWS+'_Inty',0,True)
575                if (Plot == 'FwHm' or Plot == 'All'):
576                        w_plot=mp.plotSpectrum(inputWS+'_FwHm',0,True)
577                if (Plot == 'Beta' or Plot == 'All'):
578                        s_plot=mp.plotSpectrum(inputWS+'_Beta',0,True)
579        if (Plot == 'Fit' or Plot == 'All'):
580                fWS = inputWS+'_Result_0'
581                f_plot=mp.plotSpectrum(fWS,res_plot,False)
582
583# Quest programs
584
585def CheckBetSig(nbs):
586        Nsig = int(nbs[1])
587        if Nsig == 0:
588                error = 'Number of sigma points is Zero'                       
589                logger.notice('ERROR *** ' + error)
590                sys.exit(error)
591        if Nsig > 200:
592                error = 'Max number of sigma points is 200'                     
593                logger.notice('ERROR *** ' + error)
594                sys.exit(error)
595        Nbet = int(nbs[0])
596        if Nbet == 0:
597                error = 'Number of beta points is Zero'                 
598                logger.notice('ERROR *** ' + error)
599                sys.exit(error)
600        if Nbet > 200:
601                error = 'Max number of beta points is 200'                     
602                logger.notice('ERROR *** ' + error)
603                sys.exit(error)
604        return Nbet,Nsig
605
606def QuestRun(samWS,resWS,nbs,erange,nbins,fitOp,Loop,Verbose,Plot,Save):
607        StartTime('Quest')
608        workdir = config['defaultsave.directory']
609        array_len = 4096                           # length of array in Fortran
610        CheckXrange(erange,'Energy')
611        nbin,nrbin = CheckBinning(nbins)
612        if Verbose:
613                logger.notice('Sample is ' + samWS)
614                logger.notice('Resolution is ' + resWS)
615        CheckAnalysers(samWS,resWS,Verbose)
616        nsam,ntc = CheckHistZero(samWS)
617        if Loop != True:
618                nsam = 1
619        efix = getEfixed(samWS)
620        theta,Q = GetThetaQ(samWS)
621        nres,ntr = CheckHistZero(resWS)
622        if nres == 1:
623                prog = 'Qst'                        # res file
624        else:
625                error = 'Stretched Exp ONLY works with RES file'                       
626                logger.notice('ERROR *** ' + error)
627                sys.exit(error)
628        if Verbose:
629                logger.notice(' Number of spectra = '+str(nsam))
630                logger.notice(' Erange : '+str(erange[0])+' to '+str(erange[1]))
631        dtn,xsc = ReadNormFile(fitOp[3],nsam,Verbose)
632        fname = samWS[:-4] + '_'+ prog
633        wrks=workdir + samWS[:-4]
634        if Verbose:
635                logger.notice(' lptfile : ' + wrks +'_Qst.lpt')
636        lwrk=len(wrks)
637        wrks.ljust(140,' ')
638        wrkr=resWS
639        wrkr.ljust(140,' ')
640        wrk = [wrks, wrkr]
641        Nbet,Nsig = CheckBetSig(nbs)
642        eBet0 = np.zeros(Nbet)                  # set errors to zero
643        eSig0 = np.zeros(Nsig)                  # set errors to zero
644        rscl = 1.0
645        Qaxis = ''
646        for m in range(0,nsam):
647                if Verbose:
648                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
649                nsp = m+1
650                nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
651                Ndat = nout[0]
652                Imin = nout[1]
653                Imax = nout[2]
654                Nb,Xb,Yb,Eb = GetXYE(resWS,0,array_len)
655                numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin, Nbet, Nsig]
656                reals = [efix, theta[m], rscl, bnorm]
657                xsout,ysout,xbout,ybout,zpout=Que.quest(numb,Xv,Yv,Ev,reals,fitOp,
658                                                                                        Xdat,Xb,Yb,wrks,wrkr,lwrk)
659                dataXs = xsout[:Nsig]               # reduce from fixed Fortran array
660                dataYs = ysout[:Nsig]
661                dataXb = xbout[:Nbet]
662                dataYb = ybout[:Nbet]
663                zpWS = fname + '_Zp' +str(m)
664                if (m > 0):
665                        Qaxis += ','
666                Qaxis += str(Q[m])
667                for n in range(0,Nsig):
668                        yfit_list = np.split(zpout[:Nsig*Nbet],Nsig)
669                        dataYzp = yfit_list[n]
670                        if n == 0:
671                                CreateWorkspace(OutputWorkspace=zpWS, DataX=xbout[:Nbet], DataY=dataYzp[:Nbet], DataE=eBet0,
672                                        Nspec=1, UnitX='MomentumTransfer')
673                        else:
674                                CreateWorkspace(OutputWorkspace='__Zpt', DataX=xbout[:Nbet], DataY=dataYzp[:Nbet], DataE=eBet0,
675                                        Nspec=1, UnitX='MomentumTransfer')
676                                ConjoinWorkspaces(InputWorkspace1=zpWS, InputWorkspace2='__Zpt', CheckOverlapping=False)                               
677                if m == 0:
678                        xSig = dataXs
679                        ySig = dataYs
680                        eSig = eSig0           
681                        xBet = dataXb
682                        yBet = dataYb
683                        eBet = eBet0           
684                        groupZ = zpWS
685                else:
686                        xSig = np.append(xSig,dataXs)
687                        ySig = np.append(ySig,dataYs)
688                        eSig = np.append(eSig,eSig0)
689                        xBet = np.append(xBet,dataXb)
690                        yBet = np.append(yBet,dataYb)
691                        eBet = np.append(eBet,eBet0)
692                        groupZ = groupZ +','+ zpWS
693        CreateWorkspace(OutputWorkspace=fname+'_Sigma', DataX=xSig, DataY=ySig, DataE=eSig,
694                Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
695        CreateWorkspace(OutputWorkspace=fname+'_Beta', DataX=xBet, DataY=yBet, DataE=eBet,
696                Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
697        group = fname + '_Sigma,'+ fname + '_Beta'
698        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_Fit')     
699        if Loop:
700                GroupWorkspaces(InputWorkspaces=groupZ,OutputWorkspace=fname+'_Contour')
701        if Save:
702                fpath = os.path.join(workdir,fname+'_Fit.nxs')
703                SaveNexusProcessed(InputWorkspace=fname+'_Fit', Filename=fpath)
704                cpath = os.path.join(workdir,fname+'_Contour.nxs')
705                SaveNexusProcessed(InputWorkspace=fname+'_Contour', Filename=cpath)
706                if Verbose:
707                        logger.notice('Output file for Fit : ' + fpath)
708                        logger.notice('Output file for Contours : ' + cpath)
709        if (Plot != 'None'):
710                QuestPlot(fname,Plot)
711        EndTime('Quest')
712
713def QuestPlot(inputWS,Plot):
714        if (Plot == 'Sigma' or Plot == 'All'):
715                s_graph = mp.importMatrixWorkspace(inputWS+'_Sigma').plotGraph2D()
716        s_layer = s_graph.activeLayer().setAxisTitle(2, 'Sigma')
717        if (Plot == 'Beta' or Plot == 'All'):
718                b_graph = mp.importMatrixWorkspace(inputWS+'_Beta').plotGraph2D()
719        b_layer = b_graph.activeLayer().setAxisTitle(2, 'Beta')
720
721# ResNorm programs
722
723def ResNormRun(vname,rname,erange,nbins,Verbose,Plot,Save):
724        StartTime('ResNorm')
725        workdir = config['defaultsave.directory']
726        array_len = 4096                                    # length of Fortran array
727        CheckXrange(erange,'Energy')
728        nbin,nrbin = CheckBinning(nbins)
729        CheckAnalysers(vname,rname,Verbose)
730        nvan,ntc = CheckHistZero(vname)
731        theta,Q = GetThetaQ(vname)
732        efix = getEfixed(vname)
733        nres,ntr = CheckHistZero(rname)
734        nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(vname,0,erange,nbin)
735        Ndat = nout[0]
736        Imin = nout[1]
737        Imax = nout[2]
738        wrks=workdir + vname[:-4]
739        if Verbose:
740                logger.notice(' Number of spectra = '+str(nvan))
741                logger.notice(' lptfile : ' + wrks +'_resnrm.lpt')
742        lwrk=len(wrks)
743        wrks.ljust(140,' ')                              # pad for fioxed Fortran length
744        wrkr=rname
745        wrkr.ljust(140,' ')
746        Nb,Xb,Yb,Eb = GetXYE(rname,0,array_len)
747        rscl = 1.0
748        xPar = np.array([theta[0]])
749        for m in range(1,nvan):
750                xPar = np.append(xPar,theta[m])
751        ePar = np.zeros(nvan)
752        fname = vname[:-4]
753        for m in range(0,nvan):
754                if Verbose:
755                        logger.notice('Group ' +str(m)+ ' at angle '+ str(theta[m]))
756                ntc,Xv,Yv,Ev = GetXYE(vname,m,array_len)
757                nsp = m+1
758                numb = [nvan, nsp, ntc, Ndat, nbin, Imin, Imax, Nb]
759                reals = [efix, theta[0], rscl, bnorm]
760                nd,xout,yout,eout,yfit,pfit=resnorm.resnorm(numb,Xv,Yv,Ev,reals,
761                                                                        Xdat,Xb,Yb,wrks,wrkr,lwrk)
762                if Verbose:
763                        message = ' Fit paras : '+str(pfit[0])+' '+str(pfit[1])
764                        logger.notice(message)
765                dataX = xout[:nd]
766                dataX = np.append(dataX,2*xout[nd-1]-xout[nd-2])
767                if m == 0:
768                        yPar1 = np.array([pfit[0]])
769                        yPar2 = np.array([pfit[1]])
770                        CreateWorkspace(OutputWorkspace='Data', DataX=dataX, DataY=yout[:nd], DataE=eout[:nd],
771                                NSpec=1, UnitX='DeltaE')
772                        CreateWorkspace(OutputWorkspace='Fit', DataX=dataX, DataY=yfit[:nd], DataE=np.zeros(nd),
773                                NSpec=1, UnitX='DeltaE')
774                else:
775                        yPar1 = np.append(yPar1,pfit[0])
776                        yPar2 = np.append(yPar2,pfit[1])
777                        CreateWorkspace(OutputWorkspace='__datmp', DataX=dataX, DataY=yout[:nd], DataE=eout[:nd],
778                                NSpec=1, UnitX='DeltaE')
779                        ConjoinWorkspaces(InputWorkspace1='Data', InputWorkspace2='__datmp', CheckOverlapping=False)                           
780                        CreateWorkspace(OutputWorkspace='__f1tmp', DataX=dataX, DataY=yfit[:nd], DataE=np.zeros(nd),
781                                NSpec=1, UnitX='DeltaE')
782                        ConjoinWorkspaces(InputWorkspace1='Fit', InputWorkspace2='__f1tmp', CheckOverlapping=False)                             
783        CreateWorkspace(OutputWorkspace=fname+'_ResNorm_Intensity', DataX=xPar, DataY=yPar1, DataE=xPar,
784                NSpec=1, UnitX='MomentumTransfer')
785        CreateWorkspace(OutputWorkspace=fname+'_ResNorm_Stretch', DataX=xPar, DataY=yPar2, DataE=xPar,
786                NSpec=1, UnitX='MomentumTransfer')
787        group = fname + '_ResNorm_Intensity,'+ fname + '_ResNorm_Stretch'
788        GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname+'_ResNorm_Paras')
789        GroupWorkspaces(InputWorkspaces='Data,Fit',OutputWorkspace=fname+'_ResNorm_Fit')
790        if Save:
791                par_path = os.path.join(workdir,fname+'_ResNorm_Paras.nxs')
792                SaveNexusProcessed(InputWorkspace=fname+'_ResNorm_Paras', Filename=par_path)
793                fit_path = os.path.join(workdir,fname+'_ResNorm_Fit.nxs')
794                SaveNexusProcessed(InputWorkspace=fname+'_ResNorm_Fit', Filename=fit_path)
795                if Verbose:
796                        logger.notice('Parameter file created : ' + par_path)
797                        logger.notice('Fit file created : ' + fit_path)
798        if (Plot != 'None'):
799                ResNormPlot(fname,Plot)
800        EndTime('ResNorm')
801
802def ResNormPlot(inputWS,Plot):
803        if (Plot == 'Intensity' or Plot == 'All'):
804                iWS = inputWS + '_ResNorm_Intensity'
805                i_plot=mp.plotSpectrum(iWS,0,False)
806        if (Plot == 'Stretch' or Plot == 'All'):
807                sWS = inputWS + '_ResNorm_Stretch'
808                s_plot=mp.plotSpectrum(sWS,0,False)
809        if (Plot == 'Fit' or Plot == 'All'):
810                fWS = inputWS + '_ResNorm_Fit'
811                f_plot=mp.plotSpectrum(fWS,0,False)
812
813# Jump programs
814
815def JumpRun(samWS,jump,Verbose,Plot,Save):
816# Chudley-Elliott    HWHM=A*(1-sin*(Q*K)/(Q*K))
817# for Q->0 W=A*Q^2*K^2/6
818# Singwi-Sjolander  HWHM=A*(1-exp(-r*Q^2))
819# for Q->0 W=A*Q^2*r
820        StartTime('Jump fit : '+jump+' ; ')
821        workdir = config['defaultsave.directory']
822        if Verbose:
823                logger.notice('Parameters in ' + samWS+'_Workspace')
824        x = mtd[samWS].readX(0)
825        xmax = x[len(x)-1]
826        if jump == 'CE':
827                aval = xmax
828                bval = 1.5
829                func = 'name=UserFunction, Formula=a*(1-(sin(x*b))/(x*b)), a='+str(aval)+', b='+str(bval)
830        if jump == 'SS':
831                aval = xmax
832                bval = 0.24
833                func = 'name=UserFunction, Formula=a*(1-exp(-x*x*b)), a='+str(aval)+', b='+str(bval)
834        fitWS = samWS +'_'+jump +'fit'
835        Fit(Function=func, InputWorkspace=samWS, CreateOutput=True, Output=fitWS)
836        if Save:
837                fit_path = os.path.join(workdir,fitWS+'.nxs')
838                SaveNexusProcessed(InputWorkspace=fitWS, Filename=fit_path)
839                if Verbose:
840                        logger.notice('Fit file is ' + fit_path)
841        if Plot:
842                JumpPlot(fitWS+'_Workspace')
843        EndTime('Jump fit : '+jump+' ; ')
844
845def JumpPlot(inputWS):
846    j_plot=mp.plotSpectrum(inputWS,[0],True)
847    mp.mergePlots(j_plot,mp.plotSpectrum(inputWS,[1,2],False))