Ticket #7965: IndirectBayes.py

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

Updated file: Lines 296 & 300 modified, 303-305 added.

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