Context Navigation

Back to Ticket #8585

Ticket #8585: IndirectBayes.py

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

Download in other formats:

Original Format