Context Navigation

Back to Ticket #6426

Ticket #6426: IndirectDataAnalysis.py

File IndirectDataAnalysis.py, 35.3 KB (added by Martyn Gigg, 8 years ago)

Line
1	from mantid.simpleapi import *
2	from IndirectImport import import_mantidplot
3	mp = import_mantidplot()
4	from IndirectCommon import *
5	from mantid import config, logger
6	import math, re, os.path, numpy as np
7
8	##############################################################################
9	# Misc. Helper Functions
10	##############################################################################
11
12	def concatWSs(workspaces, unit, name):
13	dataX = []
14	dataY = []
15	dataE = []
16	for ws in workspaces:
17	readX = mtd[ws].readX(0)
18	readY = mtd[ws].readY(0)
19	readE = mtd[ws].readE(0)
20	for i in range(0, len(readX)):
21	dataX.append(readX[i])
22	for i in range(0, len(readY)):
23	dataY.append(readY[i])
24	dataE.append(readE[i])
25	CreateWorkspace(OutputWorkspace=name, DataX=dataX, DataY=dataY, DataE=dataE,
26	NSpec=len(workspaces), UnitX=unit)
27
28	def split(l, n):
29	#Yield successive n-sized chunks from l.
30	for i in xrange(0, len(l), n):
31	yield l[i:i+n]
32
33	def segment(l, fromIndex, toIndex):
34	for i in xrange(fromIndex, toIndex + 1):
35	yield l[i]
36
37	def trimData(nSpec, vals, min, max):
38	result = []
39	chunkSize = len(vals) / nSpec
40	assert min >= 0, 'trimData: min is less then zero'
41	assert max <= chunkSize - 1, 'trimData: max is greater than the number of spectra'
42	assert min <= max, 'trimData: min is greater than max'
43	chunks = split(vals,chunkSize)
44	for chunk in chunks:
45	seg = segment(chunk,min,max)
46	for val in seg:
47	result.append(val)
48	return result
49
50	##############################################################################
51	# ConvFit
52	##############################################################################
53
54	def confitParsToWS(Table, Data, BackG='FixF', specMin=0, specMax=-1):
55	if ( specMax == -1 ):
56	specMax = mtd[Data].getNumberHistograms() - 1
57	dataX = createQaxis(Data)
58	xAxisVals = []
59	xAxisTrimmed = []
60	dataY = []
61	dataE = []
62	names = ''
63	ws = mtd[Table]
64	cName = ws.getColumnNames()
65	nSpec = ( ws.columnCount() - 1 ) / 2
66	for spec in range(0,nSpec):
67	yCol = (spec*2)+1
68	yAxis = cName[(spec*2)+1]
69	if re.search('HWHM$', yAxis) or re.search('Height$', yAxis):
70	xAxisVals += dataX
71	if (len(names) > 0):
72	names += ","
73	names += yAxis
74	eCol = (spec*2)+2
75	eAxis = cName[(spec*2)+2]
76	for row in range(0, ws.rowCount()):
77	dataY.append(ws.cell(row,yCol))
78	dataE.append(ws.cell(row,eCol))
79	else:
80	nSpec -= 1
81	outNm = Table + "_Workspace"
82	xAxisTrimmed = trimData(nSpec, xAxisVals, specMin, specMax)
83	CreateWorkspace(OutputWorkspace=outNm, DataX=xAxisTrimmed, DataY=dataY, DataE=dataE,
84	Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
85	VerticalAxisValues=names)
86	return outNm
87
88	def confitPlotSeq(inputWS, plot):
89	nhist = mtd[inputWS].getNumberHistograms()
90	if ( plot == 'All' ):
91	mp.plotSpectrum(inputWS, range(0, nhist), True)
92	return
93	plotSpecs = []
94	if ( plot == 'Intensity' ):
95	res = 'Height$'
96	elif ( plot == 'HWHM' ):
97	res = 'HWHM$'
98	for i in range(0,nhist):
99	title = mtd[inputWS].getAxis(1).label(i)
100	if re.search(res, title):
101	plotSpecs.append(i)
102	mp.plotSpectrum(inputWS, plotSpecs, True)
103
104	def confitSeq(inputWS, func, startX, endX, save, plot, ftype, bg, specMin, specMax, Verbose=True):
105	StartTime('ConvFit')
106	workdir = config['defaultsave.directory']
107	input = inputWS+',i' + str(specMin)
108	if (specMax == -1):
109	specMax = mtd[inputWS].getNumberHistograms() - 1
110	for i in range(specMin + 1, specMax + 1):
111	input += ';'+inputWS+',i'+str(i)
112	(instr, run) = getInstrRun(inputWS)
113	run_name = instr + run
114	outNm = getWSprefix(inputWS) + 'conv_' + ftype + bg + str(specMin) + "_to_" + str(specMax)
115	if Verbose:
116	logger.notice(func)
117	PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func,
118	StartX=startX, EndX=endX, FitType='Sequential')
119	wsname = confitParsToWS(outNm, inputWS, bg, specMin, specMax)
120	RenameWorkspace(InputWorkspace=outNm,
121	OutputWorkspace=outNm + "_Parameters")
122	if save:
123	SaveNexusProcessed(InputWorkspace=wsname, Filename=wsname+'.nxs')
124	if plot != 'None':
125	confitPlotSeq(wsname, plot)
126	EndTime('ConvFit')
127
128	##############################################################################
129	# Elwin
130	##############################################################################
131
132	def elwin(inputFiles, eRange, Save=False, Verbose=True, Plot=False):
133	StartTime('ElWin')
134	Verbose = True
135	workdir = config['defaultsave.directory']
136	CheckXrange(eRange,'Energy')
137	tempWS = '__temp'
138	if Verbose:
139	range1 = str(eRange[0])+' to '+str(eRange[1])
140	if ( len(eRange) == 4 ):
141	range2 = str(eRange[2])+' to '+str(eRange[3])
142	logger.notice('Using 2 energy ranges from '+range1+' & '+range2)
143	elif ( len(eRange) == 2 ):
144	logger.notice('Using 1 energy range from '+range1)
145	nr = 0
146	inputRuns = sorted(inputFiles)
147	Vrun = []
148	for file in inputRuns:
149	(direct, file_name) = os.path.split(file)
150	(root, ext) = os.path.splitext(file_name)
151	LoadNexus(Filename=file, OutputWorkspace=tempWS)
152	if Verbose:
153	logger.notice('Reading file : '+file)
154	nsam,ntc = CheckHistZero(tempWS)
155	if ( len(eRange) == 4 ):
156	ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1],
157	Range2Start=eRange[2], Range2End=eRange[3],
158	OutputInQ='__eq1', OutputInQSquared='__eq2')
159	elif ( len(eRange) == 2 ):
160	ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1],
161	OutputInQ='__eq1', OutputInQSquared='__eq2')
162	(instr, last) = getInstrRun(root)
163	q1 = np.array(mtd['__eq1'].readX(0))
164	i1 = np.array(mtd['__eq1'].readY(0))
165	e1 = np.array(mtd['__eq1'].readE(0))
166	q2 = np.array(mtd['__eq2'].readX(0))
167	inY = mtd['__eq2'].readY(0)
168	inE = mtd['__eq2'].readE(0)
169	logy = []
170	loge = []
171	for i in range(0, len(inY)):
172	if(inY[i] == 0):
173	ly = math.log(0.000000000001)
174	else:
175	ly = math.log(inY[i])
176	logy.append(ly)
177	if( inY[i]+inE[i] == 0 ):
178	le = math.log(0.000000000001)-ly
179	else:
180	le = math.log(inY[i]+inE[i])-ly
181	loge.append(le)
182	i2 = np.array(logy)
183	e2 = np.array(loge)
184	Vrun.append(float(last))
185	if (nr == 0):
186	first = getWSprefix(tempWS,root)
187	datX1 = q1
188	datY1 = i1
189	datE1 = e1
190	datX2 = q2
191	datY2 = i2
192	datE2 = e2
193	Vaxis = last
194	else:
195	datX1 = np.append(datX1,q1)
196	datY1 = np.append(datY1,i1)
197	datE1 = np.append(datE1,e1)
198	datX2 = np.append(datX2,q2)
199	datY2 = np.append(datY2,i2)
200	datE2 = np.append(datE2,e2)
201	Vaxis += ','+last
202	nr += 1
203	DeleteWorkspace(tempWS)
204	DeleteWorkspace('__eq1')
205	DeleteWorkspace('__eq2')
206	if (nr == 1):
207	ename = first[:-1]
208	else:
209	ename = first+'to_'+last
210	e1WS = ename+'_eq1'
211	e2WS = ename+'_eq2'
212	elwWS = ename+'_elw' # temporary fix to do plotting
213	CreateWorkspace(OutputWorkspace=elwWS, DataX=datX1, DataY=datY1, DataE=datE1,
214	Nspec=nr, UnitX='MomentumTransfer')
215	CreateWorkspace(OutputWorkspace=e1WS, DataX=datX1, DataY=datY1, DataE=datE1,
216	Nspec=nr, UnitX='MomentumTransfer', VerticalAxisUnit='Text', VerticalAxisValues=Vaxis)
217	replace_workspace_axis(e1WS, Vrun, '')
218	CreateWorkspace(OutputWorkspace=e2WS, DataX=datX2, DataY=datY2, DataE=datE2,
219	Nspec=nr, UnitX='QSquared', VerticalAxisUnit='Text', VerticalAxisValues=Vaxis)
220	replace_workspace_axis(e2WS, Vrun, '')
221	if Save:
222	e1_path = os.path.join(workdir, e1WS+'.nxs') # path name for nxs file
223	e2_path = os.path.join(workdir, e2WS+'.nxs') # path name for nxs file
224	if Verbose:
225	logger.notice('Creating file : '+e1_path)
226	logger.notice('Creating file : '+e2_path)
227	SaveNexusProcessed(InputWorkspace=e1WS, Filename=e1_path)
228	SaveNexusProcessed(InputWorkspace=e2WS, Filename=e2_path)
229	if Plot:
230	elwinPlot(e1WS,e2WS)
231	EndTime('Elwin')
232	return e1WS,e2WS
233
234	def elwinPlot(eq1,eq2):
235	nhist = mtd[eq1].getNumberHistograms() # no. of hist/groups in sam
236	nBins = mtd[eq1].blocksize()
237	lastXeq1 = mtd[eq1].readX(0)[nBins-1]
238	graph1 = mp.plotSpectrum(eq1, range(0,nhist))
239	layer1 = graph1.activeLayer()
240	layer1.setScale(mp.Layer.Bottom, 0.0, lastXeq1)
241	layer1.setAxisTitle(mp.Layer.Left,'Elastic Intensity')
242	nBins = mtd[eq2].blocksize()
243	lastXeq2 = mtd[eq2].readX(0)[nBins-1]
244	graph2 = mp.plotSpectrum(eq2, range(0,nhist))
245	layer2 = graph2.activeLayer()
246	layer2.setScale(mp.Layer.Bottom, 0.0, lastXeq2)
247	layer2.setAxisTitle(mp.Layer.Left,'log(Elastic Intensity)')
248
249	##############################################################################
250	# Fury
251	##############################################################################
252
253	def furyPlot(inWS, spec):
254	graph = mp.plotSpectrum(inWS, spec)
255	layer = graph.activeLayer()
256	layer.setScale(mp.Layer.Left, 0, 1.0)
257
258	def fury(sam_files, res_file, rebinParam, RES=True, Save=False, Verbose=False,
259	Plot=False):
260	StartTime('Fury')
261	Verbose = True
262	workdir = config['defaultsave.directory']
263	LoadNexus(Filename=sam_files[0], OutputWorkspace='__sam_tmp') # SAMPLE
264	nsam,npt = CheckHistZero('__sam_tmp')
265	Xin = mtd['__sam_tmp'].readX(0)
266	d1 = Xin[1]-Xin[0]
267	if d1 < 1e-8:
268	error = 'Data energy bin is zero'
269	logger.notice('ERROR *** ' + error)
270	sys.exit(error)
271	d2 = Xin[npt-1]-Xin[npt-2]
272	dmin = min(d1,d2)
273	pars = rebinParam.split(',')
274	if (float(pars[1]) <= dmin):
275	error = 'EWidth = ' + pars[1] + ' < smallest Eincr = ' + str(dmin)
276	logger.notice('ERROR *** ' + error)
277	sys.exit(error)
278	outWSlist = []
279	# Process RES Data Only Once
280	if Verbose:
281	logger.notice('Reading RES file : '+res_file)
282	LoadNexus(Filename=res_file, OutputWorkspace='res_data') # RES
283	CheckAnalysers('__sam_tmp','res_data',Verbose)
284	nres,nptr = CheckHistZero('res_data')
285	if nres > 1:
286	CheckHistSame('__sam_tmp','Sample','res_data','Resolution')
287	DeleteWorkspace('__sam_tmp')
288	Rebin(InputWorkspace='res_data', OutputWorkspace='res_data', Params=rebinParam)
289	Integration(InputWorkspace='res_data', OutputWorkspace='res_int')
290	ConvertToPointData(InputWorkspace='res_data', OutputWorkspace='res_data')
291	ExtractFFTSpectrum(InputWorkspace='res_data', OutputWorkspace='res_fft', FFTPart=2)
292	Divide(LHSWorkspace='res_fft', RHSWorkspace='res_int', OutputWorkspace='res')
293	for sam_file in sam_files:
294	(direct, filename) = os.path.split(sam_file)
295	(root, ext) = os.path.splitext(filename)
296	if (ext == '.nxs'): # input is file
297	if Verbose:
298	logger.notice('Reading sample file : '+sam_file)
299	LoadNexus(Filename=sam_file, OutputWorkspace='sam_data') # SAMPLE
300	Rebin(InputWorkspace='sam_data', OutputWorkspace='sam_data', Params=rebinParam)
301	else: # input is workspace
302	Rebin(InputWorkspace=sam_file, OutputWorkspace='sam_data', Params=rebinParam)
303	Integration(InputWorkspace='sam_data', OutputWorkspace='sam_int')
304	ConvertToPointData(InputWorkspace='sam_data', OutputWorkspace='sam_data')
305	ExtractFFTSpectrum(InputWorkspace='sam_data', OutputWorkspace='sam_fft', FFTPart=2)
306	Divide(LHSWorkspace='sam_fft', RHSWorkspace='sam_int', OutputWorkspace='sam')
307	savefile = getWSprefix('sam_data', root) + 'iqt' # Create save file name
308	outWSlist.append(savefile)
309	Divide(LHSWorkspace='sam', RHSWorkspace='res', OutputWorkspace=savefile)
310	# Cleanup Sample Files
311	DeleteWorkspace('sam_data')
312	DeleteWorkspace('sam_int')
313	DeleteWorkspace('sam_fft')
314	DeleteWorkspace('sam')
315	# Crop nonsense values off workspace
316	bin = int(math.ceil(mtd[savefile].blocksize()/2.0))
317	binV = mtd[savefile].dataX(0)[bin]
318	CropWorkspace(InputWorkspace=savefile, OutputWorkspace=savefile, XMax=binV)
319	if Save:
320	opath = os.path.join(workdir, savefile+'.nxs') # path name for nxs file
321	SaveNexusProcessed(InputWorkspace=savefile, Filename=opath)
322	if Verbose:
323	logger.notice('Output file : '+opath)
324	# Clean Up RES files
325	DeleteWorkspace('res_data')
326	DeleteWorkspace('res_int')
327	DeleteWorkspace('res_fft')
328	DeleteWorkspace('res')
329	if Plot:
330	specrange = range(0,mtd[outWSlist[0]].getNumberHistograms())
331	furyPlot(outWSlist, specrange)
332	EndTime('Fury')
333	return outWSlist
334
335	##############################################################################
336	# FuryFit
337	##############################################################################
338
339	def furyfitParsToWS(Table, Data, option):
340	nopt = len(option)
341	if nopt == 2:
342	npeak = option[0]
343	type = option[1]
344	elif nopt == 4:
345	npeak = 2
346	type = 'SE'
347	else:
348	logger.notice('Bad option : ' +option)
349	Q = createQaxis(Data)
350	nQ = len(Q)
351	ws = mtd[Table]
352	rCount = ws.rowCount()
353	cCount = ws.columnCount()
354	cName = ws.getColumnNames()
355	Qa = np.array(Q)
356	A0v = ws.column(1) #bgd value
357	A0e = ws.column(2) #bgd error
358	Iy1 = ws.column(5) #intensity1 value
359	Ie1 = ws.column(2) #intensity1 error = bgd
360	dataX = Qa
361	dataY = np.array(A0v)
362	dataE = np.array(A0e)
363	names = cName[1]
364	dataX = np.append(dataX,Qa)
365	dataY = np.append(dataY,np.array(Iy1))
366	dataE = np.append(dataE,np.array(Ie1))
367	names += ","+cName[5]
368	Ty1 = ws.column(7) #tau1 value
369	Te1 = ws.column(8) #tau1 error
370	dataX = np.append(dataX,Qa)
371	dataY = np.append(dataY,np.array(Ty1))
372	dataE = np.append(dataE,np.array(Te1))
373	names += ","+cName[7]
374	nSpec = 3
375	if npeak == 1:
376	By1 = ws.column(9) #beta1 value
377	Be1 = ws.column(10) #beta2 error
378	dataX = np.append(dataX,Qa)
379	dataY = np.append(dataY,np.array(By1))
380	dataE = np.append(dataE,np.array(Be1))
381	names += ","+cName[9]
382	nSpec += 1
383	if npeak == 2:
384	Iy2 = ws.column(9) #intensity2 value
385	Ie2 = ws.column(10) #intensity2 error
386	dataX = np.append(dataX,Qa)
387	dataY = np.append(dataY,np.array(Iy2))
388	dataE = np.append(dataE,np.array(Ie2))
389	names += ","+cName[9]
390	nSpec += 1
391	Ty2 = ws.column(11) #tau2 value
392	Te2 = ws.column(12) #tau2 error
393	dataX = np.append(dataX,Qa)
394	dataY = np.append(dataY,np.array(Ty2))
395	dataE = np.append(dataE,np.array(Te2))
396	names += ","+cName[11]
397	nSpec += 1
398	wsname = Table + "_Workspace"
399	CreateWorkspace(OutputWorkspace=wsname, DataX=dataX, DataY=dataY, DataE=dataE,
400	Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
401	VerticalAxisValues=names)
402	return wsname
403
404	def furyfitPlotSeq(inputWS, Plot):
405	nHist = mtd[inputWS].getNumberHistograms()
406	if ( Plot == 'All' ):
407	mp.plotSpectrum(inputWS, range(0, nHist), True)
408	return
409	plotSpecs = []
410	if ( Plot == 'Intensity' ):
411	res = 'Intensity$'
412	if ( Plot == 'Tau' ):
413	res = 'Tau$'
414	elif ( Plot == 'Beta' ):
415	res = 'Beta$'
416	for i in range(0, nHist):
417	title = mtd[inputWS].getAxis(1).label(i)
418	if ( re.search(res, title) ):
419	plotSpecs.append(i)
420	mp.plotSpectrum(inputWS, plotSpecs, True)
421
422	def furyfitSeq(inputWS, func, ftype, startx, endx, Save, Plot, Verbose = True):
423	StartTime('FuryFit')
424	workdir = config['defaultsave.directory']
425	input = inputWS+',i0'
426	nHist = mtd[inputWS].getNumberHistograms()
427	for i in range(1,nHist):
428	input += ';'+inputWS+',i'+str(i)
429	outNm = getWSprefix(inputWS) + 'fury_' + ftype + "0_to_" + str(nHist-1)
430	option = ftype[:-2]
431	if Verbose:
432	logger.notice('Option: '+option)
433	logger.notice(func)
434	PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func,
435	StartX=startx, EndX=endx, FitType='Sequential')
436	wsname = furyfitParsToWS(outNm, inputWS, option)
437	RenameWorkspace(InputWorkspace=outNm, OutputWorkspace=outNm+"_Parameters")
438	if Save:
439	opath = os.path.join(workdir, wsname+'.nxs') # path name for nxs file
440	SaveNexusProcessed(InputWorkspace=wsname, Filename=opath)
441	if Verbose:
442	logger.notice('Output file : '+opath)
443	if ( Plot != 'None' ):
444	furyfitPlotSeq(wsname, Plot)
445	EndTime('FuryFit')
446	return mtd[wsname]
447
448	def furyfitMultParsToWS(Table, Data):
449	dataX = []
450	dataA0v = []
451	dataA0e = []
452	dataY1 = []
453	dataE1 = []
454	dataY2 = []
455	dataE2 = []
456	dataY3 = []
457	dataE3 = []
458	ws = mtd[Table+'_Parameters']
459	rCount = ws.rowCount()
460	cCount = ws.columnCount()
461	logger.notice(' Cols : '+str(cCount))
462	nSpec = ( rCount - 1 ) / 5
463	for spec in range(0,nSpec):
464	A0val = 1
465	A0err = 2
466	ival = 5 #intensity value
467	ierr = 6 #intensity error
468	tval = 7 #tau value
469	terr = 8 #tau error
470	bval = 9 #beta value
471	bval = 10 #beta error
472	dataX.append(spec)
473	dataA0v.append(ws.cell(spec,A0val))
474	dataA0e.append(ws.cell(spec,A0err))
475	dataY1.append(ws.cell(spec,ival))
476	dataE1.append(ws.cell(spec,A0err))
477	dataY2.append(ws.cell(spec,tval))
478	dataE2.append(ws.cell(spec,terr))
479	dataY3.append(ws.cell(spec,bval))
480	dataE3.append(ws.cell(spec,berr))
481	suffix = 'S'
482	wsname = Table + '_' + suffix
483	CreateWorkspace(OutputWorkspace=wsname, DataX=dataX, DataY=dataY1, DataE=dataE1,
484	Nspec=1, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
485	VerticalAxisValues='Intensity')
486	CreateWorkspace(OutputWorkspace='__multmp', DataX=dataX, DataY=dataY2, DataE=dataE2,
487	Nspec=1, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
488	VerticalAxisValues='Tau')
489	ConjoinWorkspaces(InputWorkspace1=wsname, InputWorkspace2='__multmp', CheckOverlapping=False)
490	CreateWorkspace(OutputWorkspace='__multmp', DataX=dataX, DataY=dataY3, DataE=dataE3,
491	Nspec=1, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
492	VerticalAxisValues='Beta')
493	ConjoinWorkspaces(InputWorkspace1=wsname, InputWorkspace2='__multmp', CheckOverlapping=False)
494	return wsname
495
496	def furyfitPlotMult(inputWS, Plot):
497	nHist = mtd[inputWS].getNumberHistograms()
498	if ( Plot == 'All' ):
499	mp.plotSpectrum(inputWS, range(0, nHist))
500	return
501	plotSpecs = []
502	if ( Plot == 'Intensity' ):
503	mp.plotSpectrum(inputWS, 0, True)
504	if ( Plot == 'Tau' ):
505	mp.plotSpectrum(inputWS, 1, True)
506	elif ( Plot == 'Beta' ):
507	mp.plotSpectrum(inputWS, 2, True)
508
509	def furyfitMult(inputWS, func, startx, endx, Save, Plot):
510	StartTime('FuryFit Mult')
511	Verbose = True
512	workdir = config['defaultsave.directory']
513	input = inputWS+',i0'
514	nHist = mtd[inputWS].getNumberHistograms()
515	for i in range(1,nHist):
516	input += ';'+inputWS+',i'+str(i)
517	outNm = getWSprefix(inputWS) + 'fury'
518	f1 = """(
519	composite=CompositeFunctionMW,Workspace=$WORKSPACE$,WSParam=(WorkspaceIndex=$INDEX$);
520	name=LinearBackground,A0=0,A1=0,ties=(A1=0);
521	name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity=1.0,Tau=0.1,Beta=1;ties=(f1.Intensity=1-f0.A0)
522	);
523	""".replace('$WORKSPACE$',inputWS)
524	func= 'composite=MultiBG;'
525	ties='ties=('
526	for i in range(0,nHist):
527	func+=f1.replace('$INDEX$',str(i))
528	if i > 0:
529	ties += 'f' + str(i) + '.f1.Beta=f0.f1.Beta'
530	if i < nHist-1:
531	ties += ','
532	ties+=')'
533	func += ties
534	logger.notice(func)
535	Fit(InputWorkspace=inputWS,Function=func,Output=outNm)
536	wsname = furyfitMultParsToWS(outNm, inputWS)
537	if Save:
538	opath = os.path.join(workdir, wsname+'.nxs') # path name for nxs file
539	SaveNexusProcessed(InputWorkspace=wsname, Filename=opath)
540	if Verbose:
541	logger.notice('Output file : '+opath)
542	if ( Plot != 'None' ):
543	furyfitPlotMult(wsname, Plot)
544	EndTime('FuryFit')
545
546	##############################################################################
547	# MSDFit
548	##############################################################################
549
550	def msdfitParsToWS(Table, xData):
551	dataX = xData
552	ws = mtd[Table+'_Table']
553	rCount = ws.rowCount()
554	yA0 = ws.column(1)
555	eA0 = ws.column(2)
556	yA1 = ws.column(3)
557	dataY1 = map(lambda x : -x, yA1)
558	eA1 = ws.column(4)
559	wsname = Table
560	CreateWorkspace(OutputWorkspace=wsname+'_a0', DataX=dataX, DataY=yA0, DataE=eA0,
561	Nspec=1, UnitX='')
562	CreateWorkspace(OutputWorkspace=wsname+'_a1', DataX=dataX, DataY=dataY1, DataE=eA1,
563	Nspec=1, UnitX='')
564	group = wsname+'_a0,'+wsname+'_a1'
565	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=wsname)
566	return wsname
567
568	def msdfitPlotSeq(inputWS, xlabel):
569	msd_plot = mp.plotSpectrum(inputWS+'_a1',0,True)
570	msd_layer = msd_plot.activeLayer()
571	msd_layer.setAxisTitle(mp.Layer.Bottom,xlabel)
572	msd_layer.setAxisTitle(mp.Layer.Left,'<u2>')
573
574	def msdfitPlotFits(lniWS, fitWS, n):
575	mfit_plot = mp.plotSpectrum(lniWS,n,True)
576	mfit_layer = mfit_plot.activeLayer()
577	mfit_layer.setAxisTitle(mp.Layer.Left,'log(Elastic Intensity)')
578	mp.mergePlots(mfit_plot,mp.plotSpectrum(fitWS+'_line',n,False))
579
580	def msdfit(inputs, startX, endX, Save=False, Verbose=True, Plot=True):
581	Verbose = True
582	StartTime('msdFit')
583	workdir = config['defaultsave.directory']
584	log_type = 'sample'
585	file = inputs[0]
586	(direct, filename) = os.path.split(file)
587	(root, ext) = os.path.splitext(filename)
588	(instr, first) = getInstrRun(filename)
589	if Verbose:
590	logger.notice('Reading Run : '+file)
591	LoadNexusProcessed(FileName=file, OutputWorkspace=root)
592	nHist = mtd[root].getNumberHistograms()
593	nfirst = int(first)
594	file_list = []
595	run_list = []
596	x_list = []
597	for nr in range(0, nHist):
598	nsam,ntc = CheckHistZero(root)
599	run_name = instr + str(mtd[root].getAxis(1).label(nr))
600	lnWS = run_name+'_lnI'
601	file_list.append(lnWS)
602	ExtractSingleSpectrum(InputWorkspace=root, OutputWorkspace=lnWS,
603	WorkspaceIndex=nr)
604	log_name = run_name+'_'+log_type
605	log_file = log_name+'.txt'
606	log_path = FileFinder.getFullPath(log_file)
607	if (log_path == ''):
608	logger.notice(' Run : '+run_name +' ; Temperature file not found')
609	xval = int(run_name[-3:])
610	xlabel = 'Run'
611	else:
612	logger.notice('Found '+log_path)
613	LoadLog(Workspace=root, Filename=log_path)
614	run_logs = mtd[root].getRun()
615	tmp = run_logs[log_name].value
616	temp = tmp[len(tmp)-1]
617	logger.notice(' Run : '+run_name+' ; Temperature = '+str(temp))
618	xval = temp
619	xlabel = 'Temperature (K)'
620	last = str(nr)
621	if (nr == 0):
622	first = run_name
623	run_list = lnWS
624	else:
625	run_list += ';'+lnWS
626	x_list.append(xval)
627	nr += 1
628	mname = root[:-4]
629	msdWS = mname+'_msd'
630	if Verbose:
631	logger.notice('Fitting Runs '+mname)
632	logger.notice('Q-range from '+str(startX)+' to '+str(endX))
633	function = 'name=LinearBackground, A0=0, A1=0'
634	PlotPeakByLogValue(Input=run_list, OutputWorkspace=msdWS+'_Table', Function=function,
635	StartX=startX, EndX=endX, FitType = 'Sequential')
636	msdfitParsToWS(msdWS, x_list)
637	nr = 0
638	lniWS = mname+'_lnI'
639	fitWS = mname+'_Fit'
640	a0 = mtd[msdWS+'_a0'].readY(0)
641	a1 = mtd[msdWS+'_a1'].readY(0)
642	for nr in range(0, nHist):
643	inWS = file_list[nr]
644	CropWorkspace(InputWorkspace=inWS,OutputWorkspace='__data',XMin=0.95startX,XMax=1.05endX)
645	xin = mtd['__data'].readX(0)
646	nxd = len(xin)-1
647	xd = []
648	yd = []
649	ed = []
650	for n in range(0,nxd):
651	line = a0[nr] - a1[nr]*xin[n]
652	xd.append(xin[n])
653	yd.append(line)
654	ed.append(0.0)
655	xd.append(xin[nxd])
656	CreateWorkspace(OutputWorkspace='__line', DataX=xd, DataY=yd, DataE=ed,
657	Nspec=1)
658	if (nr == 0):
659	RenameWorkspace(InputWorkspace='__data',OutputWorkspace=fitWS+'_data')
660	RenameWorkspace(InputWorkspace='__line',OutputWorkspace=fitWS+'_line')
661	else:
662	ConjoinWorkspaces(InputWorkspace1=fitWS+'_data', InputWorkspace2='__data', CheckOverlapping=False)
663	ConjoinWorkspaces(InputWorkspace1=fitWS+'_line', InputWorkspace2='__line', CheckOverlapping=False)
664	nr += 1
665	group = fitWS+'_data,'+ fitWS+'_line'
666	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS)
667	DeleteWorkspace(inWS)
668	if Plot:
669	msdfitPlotSeq(msdWS, xlabel)
670	msdfitPlotFits(root, fitWS, 0)
671	if Save:
672	msd_path = os.path.join(workdir, msdWS+'.nxs') # path name for nxs file
673	SaveNexusProcessed(InputWorkspace=msdWS, Filename=msd_path, Title=msdWS)
674	if Verbose:
675	logger.notice('Output msd file : '+msd_path)
676	EndTime('msdFit')
677	return msdWS
678
679	def plotInput(inputfiles,spectra=[]):
680	OneSpectra = False
681	if len(spectra) != 2:
682	spectra = [spectra[0], spectra[0]]
683	OneSpectra = True
684	workspaces = []
685	for file in inputfiles:
686	root = LoadNexus(Filename=file)
687	if not OneSpectra:
688	GroupDetectors(root, root,
689	DetectorList=range(spectra[0],spectra[1]+1) )
690	workspaces.append(root)
691	if len(workspaces) > 0:
692	graph = mp.plotSpectrum(workspaces,0)
693	layer = graph.activeLayer().setTitle(", ".join(workspaces))
694
695	##############################################################################
696	# Corrections
697	##############################################################################
698
699	def CubicFit(inputWS, spec, Verbose=False):
700	'''Uses the Mantid Fit Algorithm to fit a quadratic to the inputWS
701	parameter. Returns a list containing the fitted parameter values.'''
702	function = 'name=Quadratic, A0=1, A1=0, A2=0'
703	fit = Fit(Function=function, InputWorkspace=inputWS, WorkspaceIndex=spec,
704	CreateOutput=True, Output='Fit')
705	table = mtd['Fit_Parameters']
706	A0 = table.cell(0,1)
707	A1 = table.cell(1,1)
708	A2 = table.cell(2,1)
709	Abs = [A0, A1, A2]
710	if Verbose:
711	logger.notice('Group '+str(spec)+' of '+inputWS+' ; fit coefficients are : '+str(Abs))
712	return Abs
713
714	def applyCorrections(inputWS, canWS, corr, Verbose=False):
715	'''Through the PolynomialCorrection algorithm, makes corrections to the
716	input workspace based on the supplied correction values.'''
717	# Corrections are applied in Lambda (Wavelength)
718	efixed = getEfixed(inputWS) # Get efixed
719	theta,Q = GetThetaQ(inputWS)
720	sam_name = getWSprefix(inputWS)
721	ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='Wavelength',
722	EMode='Indirect', EFixed=efixed)
723	if canWS != '':
724	(instr, can_run) = getInstrRun(canWS)
725	corrections = [corr+'_1', corr+'_2', corr+'_3', corr+'_4']
726	CorrectedWS = sam_name +'Correct_'+ can_run
727	ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='Wavelength',
728	EMode='Indirect', EFixed=efixed)
729	else:
730	corrections = [corr+'_1']
731	CorrectedWS = sam_name +'Corrected'
732	nHist = mtd[inputWS].getNumberHistograms()
733	# Check that number of histograms in each corrections workspace matches
734	# that of the input (sample) workspace
735	for ws in corrections:
736	if ( mtd[ws].getNumberHistograms() != nHist ):
737	raise ValueError('Mismatch: num of spectra in '+ws+' and inputWS')
738	# Workspaces that hold intermediate results
739	CorrectedSampleWS = '__csam'
740	CorrectedCanWS = '__ccan'
741	for i in range(0, nHist): # Loop through each spectra in the inputWS
742	ExtractSingleSpectrum(InputWorkspace=inputWS, OutputWorkspace=CorrectedSampleWS,
743	WorkspaceIndex=i)
744	if ( len(corrections) == 1 ):
745	Ass = CubicFit(corrections[0], i, Verbose)
746	PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
747	Coefficients=Ass, Operation='Divide')
748	if ( i == 0 ):
749	CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
750	else:
751	ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS)
752	else:
753	ExtractSingleSpectrum(InputWorkspace=canWS, OutputWorkspace=CorrectedCanWS,
754	WorkspaceIndex=i)
755	Acc = CubicFit(corrections[3], i, Verbose)
756	PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
757	Coefficients=Acc, Operation='Divide')
758	Acsc = CubicFit(corrections[2], i, Verbose)
759	PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
760	Coefficients=Acsc, Operation='Multiply')
761	Minus(LHSWorkspace=CorrectedSampleWS, RHSWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedSampleWS)
762	Assc = CubicFit(corrections[1], i, Verbose)
763	PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
764	Coefficients=Assc, Operation='Divide')
765	if ( i == 0 ):
766	CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
767	else:
768	ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS,
769	CheckOverlapping=False)
770	ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='DeltaE',
771	EMode='Indirect', EFixed=efixed)
772	ConvertUnits(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS, Target='DeltaE',
773	EMode='Indirect', EFixed=efixed)
774	CloneWorkspace(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_rqw')
775	replace_workspace_axis(CorrectedWS+'_rqw', Q, 'MomentumTransfer')
776	RenameWorkspace(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_red')
777	if canWS != '':
778	DeleteWorkspace(CorrectedCanWS)
779	ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='DeltaE',
780	EMode='Indirect', EFixed=efixed)
781	DeleteWorkspace('Fit_NormalisedCovarianceMatrix')
782	DeleteWorkspace('Fit_Parameters')
783	DeleteWorkspace('Fit_Workspace')
784	DeleteWorkspace('corrections')
785	return CorrectedWS
786
787	def abscorFeeder(sample, container, geom, useCor):
788	'''Load up the necessary files and then passes them into the main
789	applyCorrections routine.'''
790	StartTime('ApplyCorrections')
791	Verbose = True
792	Save = True
793	PlotResult = 'Both'
794	PlotContrib = 'Spectrum'
795	workdir = config['defaultsave.directory']
796	CheckAnalysers(sample,container,Verbose)
797	s_hist,sxlen = CheckHistZero(sample)
798	sam_name = getWSprefix(sample)
799	if container != '':
800	CheckHistSame(sample,'Sample',container,'Container')
801	(instr, can_run) = getInstrRun(container)
802	if useCor:
803	if Verbose:
804	text = 'Correcting sample ' + sample
805	if container != '':
806	text += ' with ' + container
807	logger.notice(text)
808	file = sam_name + geom +'_Abs.nxs'
809	abs_path = os.path.join(workdir, file) # path name for nxs file
810	if Verbose:
811	logger.notice('Correction file :'+abs_path)
812	LoadNexus(Filename=abs_path, OutputWorkspace='corrections')
813	cor_result = applyCorrections(sample, container, 'corrections', Verbose)
814	if Save:
815	cred_path = os.path.join(workdir,cor_result+'_red.nxs')
816	SaveNexusProcessed(InputWorkspace=cor_result+'_red',Filename=cred_path)
817	if Verbose:
818	logger.notice('Output file created : '+cred_path)
819	plot_list = [cor_result+'_red',sample]
820	if ( container != '' ):
821	plot_list.append(container)
822	if (PlotResult != 'None'):
823	plotCorrResult(cor_result+'_rqw',PlotResult)
824	if (PlotContrib != 'None'):
825	plotCorrContrib(plot_list,0)
826	else:
827	if ( container == '' ):
828	sys.exit('ERROR *** Invalid options - nothing to do!')
829	else:
830	sub_result = sam_name +'Subtract_'+ can_run
831	if Verbose:
832	logger.notice('Subtracting '+container+' from '+sample)
833	Minus(LHSWorkspace=sample,RHSWorkspace=container,OutputWorkspace=sub_result)
834	CloneWorkspace(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_rqw')
835	theta,Q = GetThetaQ(sample)
836	replace_workspace_axis(sub_result+'_rqw', Q)
837	RenameWorkspace(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_red')
838	if Save:
839	sred_path = os.path.join(workdir,sub_result+'_red.nxs')
840	SaveNexusProcessed(InputWorkspace=sub_result+'_red',Filename=sred_path)
841	if Verbose:
842	logger.notice('Output file created : '+sred_path)
843	plot_list = [sub_result+'_red',sample]
844	if (Plot != 'None'):
845	plotCorrResult(sub_result+'_rqw',PlotResult)
846	if (Plot != 'None'):
847	plotCorrContrib(plot_list,0)
848	EndTime('ApplyCorrections')
849
850	def plotCorrResult(inWS,PlotResult):
851	nHist = mtd[inWS].getNumberHistograms()
852	if (Plot == 'Spectrum' or Plot == 'Both'):
853	if nHist >= 10: #only plot up to 10 hists
854	nHist = 10
855	plot_list = []
856	for i in range(0, nHist):
857	plot_list.append(i)
858	res_plot=mp.plotSpectrum(inWS,plot_list)
859	if (Plot == 'Contour' or Plot == 'Both'):
860	if nHist >= 5: #needs at least 5 hists for a contour
861	mp.importMatrixWorkspace(inWS).plotGraph2D()
862
863	def plotCorrContrib(plot_list,n):
864	con_plot=mp.plotSpectrum(plot_list,n)
865
866	def replace_workspace_axis(wsName, new_values, new_unit):
867	from mantidsimple import createNumericAxis, mtd #temporary use of old API
868	ax1 = createNumericAxis(len(new_values))
869	for i in range(len(new_values)):
870	ax1.setValue(i, new_values[i])
871	if new_unit != '':
872	ax1.setUnit(new_unit)
873	mtd[wsName].replaceAxis(1, ax1) #axis=1 is vertical

Download in other formats:

Original Format