Context Navigation

Back to Ticket #6280

Ticket #6280: IndirectBayes.py

File IndirectBayes.py, 31.1 KB (added by Nick Draper, 8 years ago)
Python routines for Bayesian analysis

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

Download in other formats:

Original Format