Context Navigation

Back to Ticket #8738

Ticket #8738: IndirectBayes.py

File IndirectBayes.py, 31.7 KB (added by Samuel Jackson, 7 years ago)
IndirectBayes, with Specner's mods

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

Download in other formats:

Original Format