Context Navigation

Back to Ticket #8388

Ticket #8388: IndirectDataAnalysis.py

File IndirectDataAnalysis.py, 50.4 KB (added by Samuel Jackson, 7 years ago)
Updated version of the script

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 split(l, n):
13	#Yield successive n-sized chunks from l.
14	for i in xrange(0, len(l), n):
15	yield l[i:i+n]
16
17	def segment(l, fromIndex, toIndex):
18	for i in xrange(fromIndex, toIndex + 1):
19	yield l[i]
20
21	def trimData(nSpec, vals, min, max):
22	result = []
23	chunkSize = len(vals) / nSpec
24	assert min >= 0, 'trimData: min is less then zero'
25	assert max <= chunkSize - 1, 'trimData: max is greater than the number of spectra'
26	assert min <= max, 'trimData: min is greater than max'
27	chunks = split(vals,chunkSize)
28	for chunk in chunks:
29	seg = segment(chunk,min,max)
30	for val in seg:
31	result.append(val)
32	return result
33
34	##############################################################################
35	# ConvFit
36	##############################################################################
37
38	def getConvFitOption(ftype, bgd, Verbose):
39	if ftype[:5] == 'Delta':
40	delta = True
41	lor = ftype[5:6]
42	else:
43	delta = False
44	lor = ftype[:1]
45	options = [bgd, delta, int(lor)]
46	if Verbose:
47	logger.notice('Fit type : Delta = ' + str(options[1]) + ' ; Lorentzians = ' + str(options[2]))
48	logger.notice('Background type : ' + options[0])
49	return options
50
51	##############################################################################
52
53	def createConvFitFun(options, par, file):
54	bgd_fun = 'name=LinearBackground,A0='
55	if options[0] == 'FixF':
56	bgd_fun = bgd_fun +str(par[0])+',A1=0,ties=(A0='+str(par[0])+',A1=0.0)'
57	if options[0] == 'FitF':
58	bgd_fun = bgd_fun +str(par[0])+',A1=0,ties=(A1=0.0)'
59	if options[0] == 'FitL':
60	bgd_fun = bgd_fun +str(par[0])+',A1='+str(par[1])
61	if options[1]:
62	ip = 3
63	else:
64	ip = 2
65	pk_1 = '(composite=Convolution;name=Resolution, FileName="'+file+'"'
66	if options[2] >= 1:
67	lor_fun = 'name=Lorentzian,Amplitude='+str(par[ip])+',PeakCentre='+str(par[ip+1])+',HWHM='+str(par[ip+2])
68	if options[2] == 2:
69	funcIndex = 1 if options[1] else 0
70	lor_2 = 'name=Lorentzian,Amplitude='+str(par[ip+3])+',PeakCentre='+str(par[ip+4])+',HWHM='+str(par[ip+5])
71	lor_fun = lor_fun +';'+ lor_2 +';ties=(f'+str(funcIndex)+'.PeakCentre=f'+str(funcIndex+1)+'.PeakCentre)'
72	if options[1]:
73	delta_fun = 'name=DeltaFunction,Height='+str(par[2])
74	lor_fun = delta_fun +';' + lor_fun
75	func = bgd_fun +';'+ pk_1 +';('+ lor_fun +'))'
76	return func
77
78	##############################################################################
79
80	def getConvFitResult(inputWS, resFile, outNm, ftype, bgd, specMin, specMax, Verbose):
81	options = getConvFitOption(ftype, bgd[:-2], Verbose)
82	params = mtd[outNm+'_Parameters']
83	A0 = params.column(1) #bgd A0 value
84	A1 = params.column(3) #bgd A1 value
85	if options[1]:
86	ip = 7
87	D1 = params.column(5) #delta value
88	else:
89	ip = 5
90	if options[2] >= 1:
91	H1 = params.column(ip) #height1 value
92	C1 = params.column(ip+2) #centre1 value
93	W1 = params.column(ip+4) #width1 value
94	if options[2] == 2:
95	H2 = params.column(ip+6) #height2 value
96	C2 = params.column(ip+8) #centre2 value
97	W2 = params.column(ip+10) #width2 value
98
99	for i in range(0,specMax-specMin):
100	paras = [A0[i], A1[i]]
101	if options[1]:
102	paras.append(D1[i])
103	if options[2] >= 1:
104	paras.append(H1[i])
105	paras.append(C1[i])
106	paras.append(W1[i])
107	if options[2] == 2:
108	paras.append(H2[i])
109	paras.append(C2[i])
110	paras.append(W2[i])
111	func = createConvFitFun(options, paras, resFile)
112	if Verbose:
113	logger.notice('Fit func : '+func)
114	fitWS = outNm + '_Result_'
115	fout = fitWS + str(i)
116	Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i+specMin,Output=fout,MaxIterations=0)
117	unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
118	unitx.setLabel('Time' , 'ns')
119	RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
120	AddSampleLog(Workspace=fout, LogName="Fit Program", LogType="String", LogText='ConvFit')
121	AddSampleLog(Workspace=fout, LogName='Background', LogType='String', LogText=str(options[0]))
122	AddSampleLog(Workspace=fout, LogName='Delta', LogType='String', LogText=str(options[1]))
123	AddSampleLog(Workspace=fout, LogName='Lorentzians', LogType='String', LogText=str(options[2]))
124	DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
125	DeleteWorkspace(fitWS+str(i)+'_Parameters')
126	if i == 0:
127	group = fout
128	else:
129	group += ',' + fout
130	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
131
132	##############################################################################
133
134	def confitParsToWS(Table, Data, specMin=0, specMax=-1):
135	if ( specMax == -1 ):
136	specMax = mtd[Data].getNumberHistograms() - 1
137	dataX = createQaxis(Data)
138	xAxisVals = []
139	xAxisTrimmed = []
140	dataY = []
141	dataE = []
142	names = ''
143	ws = mtd[Table]
144	cName = ws.getColumnNames()
145	nSpec = ( ws.columnCount() - 1 ) / 2
146	for spec in range(0,nSpec):
147	yCol = (spec*2)+1
148	yAxis = cName[(spec*2)+1]
149	if re.search('HWHM$', yAxis) or re.search('Amplitude$', yAxis):
150	xAxisVals += dataX
151	if (len(names) > 0):
152	names += ","
153	names += yAxis
154	eCol = (spec*2)+2
155	eAxis = cName[(spec*2)+2]
156	for row in range(0, ws.rowCount()):
157	dataY.append(ws.cell(row,yCol))
158	dataE.append(ws.cell(row,eCol))
159	else:
160	nSpec -= 1
161	outNm = Table + "_Workspace"
162	xAxisTrimmed = trimData(nSpec, xAxisVals, specMin, specMax)
163	CreateWorkspace(OutputWorkspace=outNm, DataX=xAxisTrimmed, DataY=dataY, DataE=dataE,
164	Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
165	VerticalAxisValues=names)
166	return outNm
167
168	##############################################################################
169
170	def confitPlotSeq(inputWS, Plot):
171	nhist = mtd[inputWS].getNumberHistograms()
172	if ( Plot == 'All' ):
173	mp.plotSpectrum(inputWS, range(0, nhist), True)
174	return
175	plotSpecs = []
176	if ( Plot == 'Intensity' ):
177	res = 'Amplitude$'
178	elif ( Plot == 'HWHM' ):
179	res = 'HWHM$'
180	for i in range(0,nhist):
181	title = mtd[inputWS].getAxis(1).label(i)
182	if re.search(res, title):
183	plotSpecs.append(i)
184	mp.plotSpectrum(inputWS, plotSpecs, True)
185
186	##############################################################################
187
188	def confitSeq(inputWS, func, startX, endX, Save, Plot, ftype, bgd, specMin, specMax, Verbose):
189	StartTime('ConvFit')
190	workdir = config['defaultsave.directory']
191	elements = func.split('"')
192	resFile = elements[1]
193	if Verbose:
194	logger.notice('Input files : '+str(inputWS))
195	input = inputWS+',i' + str(specMin)
196	if (specMax == -1):
197	specMax = mtd[inputWS].getNumberHistograms() - 1
198	for i in range(specMin + 1, specMax + 1):
199	input += ';'+inputWS+',i'+str(i)
200	(instr, run) = getInstrRun(inputWS)
201	run_name = instr + run
202	outNm = getWSprefix(inputWS) + 'conv_' + ftype + bgd + str(specMin) + "_to_" + str(specMax)
203	if Verbose:
204	logger.notice(func)
205	PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func,
206	StartX=startX, EndX=endX, FitType='Sequential')
207	wsname = confitParsToWS(outNm, inputWS, specMin, specMax)
208
209	# Add some information about convfit to the output workspace
210	options = getConvFitOption(ftype, bgd[:-2], Verbose)
211	AddSampleLog(Workspace=wsname, LogName="Fit Program", LogType="String", LogText='ConvFit')
212	AddSampleLog(Workspace=wsname, LogName='Background', LogType='String', LogText=str(options[0]))
213	AddSampleLog(Workspace=wsname, LogName='Delta', LogType='String', LogText=str(options[1]))
214	AddSampleLog(Workspace=wsname, LogName='Lorentzians', LogType='String', LogText=str(options[2]))
215
216	RenameWorkspace(InputWorkspace=outNm, OutputWorkspace=outNm + "_Parameters")
217	getConvFitResult(inputWS, resFile, outNm, ftype, bgd, specMin, specMax, Verbose)
218	if Save:
219	o_path = os.path.join(workdir, wsname+'.nxs') # path name for nxs file
220	if Verbose:
221	logger.notice('Creating file : '+o_path)
222	SaveNexusProcessed(InputWorkspace=wsname, Filename=o_path)
223	if Plot != 'None':
224	confitPlotSeq(wsname, Plot)
225	EndTime('ConvFit')
226
227	##############################################################################
228	# Elwin
229	##############################################################################
230
231	def GetTemperature(root,tempWS,log_type,Verbose):
232	(instr, run) = getInstrRun(root)
233	run_name = instr+run
234	log_name = run_name+'_'+log_type
235	run = mtd[tempWS].getRun()
236	unit1 = 'Temperature' # default values
237	unit2 = 'K'
238	if log_type in run: # test logs in WS
239	tmp = run[log_type].value
240	temp = tmp[len(tmp)-1]
241	xval = temp
242	mess = ' Run : '+run_name +' ; Temperature in log = '+str(temp)
243	else: # logs not in WS
244	logger.notice('Log parameter not found')
245	log_file = log_name+'.txt'
246	log_path = FileFinder.getFullPath(log_file)
247	if (log_path == ''): # log file does not exists
248	mess = ' Run : '+run_name +' ; Temperature file not found'
249	xval = int(run_name[-3:])
250	unit1 = 'Run-number'
251	unit2 = 'last 3 digits'
252	else: # get from log file
253	LoadLog(Workspace=tempWS, Filename=log_path)
254	run_logs = mtd[tempWS].getRun()
255	tmp = run_logs[log_name].value
256	temp = tmp[len(tmp)-1]
257	xval = temp
258	mess = ' Run : '+run_name+' ; Temperature in file = '+str(temp)
259	if Verbose:
260	logger.notice(mess)
261	unit = [unit1,unit2]
262	return xval,unit
263
264	def elwin(inputFiles, eRange, log_type='sample', Normalise = False,
265	Save=False, Verbose=False, Plot=False):
266	StartTime('ElWin')
267	workdir = config['defaultsave.directory']
268	CheckXrange(eRange,'Energy')
269	tempWS = '__temp'
270	Range2 = ( len(eRange) == 4 )
271	if Verbose:
272	range1 = str(eRange[0])+' to '+str(eRange[1])
273	if ( len(eRange) == 4 ):
274	range2 = str(eRange[2])+' to '+str(eRange[3])
275	logger.notice('Using 2 energy ranges from '+range1+' & '+range2)
276	elif ( len(eRange) == 2 ):
277	logger.notice('Using 1 energy range from '+range1)
278	nr = 0
279	inputRuns = sorted(inputFiles)
280	for file in inputRuns:
281	(direct, file_name) = os.path.split(file)
282	(root, ext) = os.path.splitext(file_name)
283	LoadNexus(Filename=file, OutputWorkspace=tempWS)
284	if (nr == 0):
285	LoadNexus(Filename=file, OutputWorkspace='run1')
286	nsam,ntc = CheckHistZero(tempWS)
287	(xval, unit) = GetTemperature(root,tempWS,log_type,Verbose)
288	if Verbose:
289	logger.notice('Reading file : '+file)
290	if ( len(eRange) == 4 ):
291	ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1],
292	Range2Start=eRange[2], Range2End=eRange[3],
293	OutputInQ='__eq1', OutputInQSquared='__eq2')
294	elif ( len(eRange) == 2 ):
295	ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1],
296	OutputInQ='__eq1', OutputInQSquared='__eq2')
297	(instr, last) = getInstrRun(root)
298	q1 = np.array(mtd['__eq1'].readX(0))
299	i1 = np.array(mtd['__eq1'].readY(0))
300	e1 = np.array(mtd['__eq1'].readE(0))
301	Logarithm(InputWorkspace='__eq2', OutputWorkspace='__eq2')
302	q2 = np.array(mtd['__eq2'].readX(0))
303	i2 = np.array(mtd['__eq2'].readY(0))
304	e2 = np.array(mtd['__eq2'].readE(0))
305	if (nr == 0):
306	CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elf')
307	first = getWSprefix(tempWS,root)
308	datX1 = q1
309	datY1 = i1
310	datE1 = e1
311	datX2 = q2
312	datY2 = i2
313	datE2 = e2
314	Tvalue = [xval]
315	Terror = [0.0]
316	Taxis = str(xval)
317	else:
318	CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elftmp')
319	ConjoinWorkspaces(InputWorkspace1='__elf', InputWorkspace2='__elftmp', CheckOverlapping=False)
320	datX1 = np.append(datX1,q1)
321	datY1 = np.append(datY1,i1)
322	datE1 = np.append(datE1,e1)
323	datX2 = np.append(datX2,q2)
324	datY2 = np.append(datY2,i2)
325	datE2 = np.append(datE2,e2)
326	Tvalue.append(xval)
327	Terror.append(0.0)
328	Taxis += ','+str(xval)
329	nr += 1
330	Txa = np.array(Tvalue)
331	Tea = np.array(Terror)
332	nQ = len(q1)
333	for nq in range(0,nQ):
334	iq = []
335	eq = []
336	for nt in range(0,len(Tvalue)):
337	ii = mtd['__elf'].readY(nt)
338	iq.append(ii[nq])
339	ie = mtd['__elf'].readE(nt)
340	eq.append(ie[nq])
341	iqa = np.array(iq)
342	eqa = np.array(eq)
343	if (nq == 0):
344	datTx = Txa
345	datTy = iqa
346	datTe = eqa
347	else:
348	datTx = np.append(datTx,Txa)
349	datTy = np.append(datTy,iqa)
350	datTe = np.append(datTe,eqa)
351	DeleteWorkspace(tempWS)
352	DeleteWorkspace('__eq1')
353	DeleteWorkspace('__eq2')
354	if (nr == 1):
355	ename = first[:-1]
356	else:
357	ename = first+'to_'+last
358
359	elfWS = ename+'_elf' # interchange Q & T
360
361	#check if temp was increasing of decreasing
362	if(datTx[0] > datTx[-1]):
363	# if so reverse data to follow natural ordering
364	datTx = datTx[::-1]
365	datTy = datTy[::-1]
366	datTe = datTe[::-1]
367
368	CreateWorkspace(OutputWorkspace=elfWS, DataX=datTx, DataY=datTy, DataE=datTe,
369	Nspec=nQ, UnitX='Energy', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q1)
370	DeleteWorkspace('__elf')
371	label = unit[0]+' / '+unit[1]
372	CopyLogs(InputWorkspace='run1', OutputWorkspace=elfWS)
373	AddSampleLog(Workspace=elfWS, LogName="vert_axis", LogType="String", LogText=label)
374	e1WS = ename+'_eq1'
375	CreateWorkspace(OutputWorkspace=e1WS, DataX=datX1, DataY=datY1, DataE=datE1,
376	Nspec=nr, UnitX='MomentumTransfer', VerticalAxisUnit='Energy', VerticalAxisValues=Taxis)
377	label = unit[0]+' / '+unit[1]
378	CopyLogs(InputWorkspace='run1', OutputWorkspace=e1WS)
379	AddSampleLog(Workspace=e1WS, LogName="vert_axis", LogType="String", LogText=label)
380	e2WS = ename+'_eq2'
381	CreateWorkspace(OutputWorkspace=e2WS, DataX=datX2, DataY=datY2, DataE=datE2,
382	Nspec=nr, UnitX='QSquared', VerticalAxisUnit='Energy', VerticalAxisValues=Taxis)
383	CopyLogs(InputWorkspace='run1', OutputWorkspace=e2WS)
384	AddSampleLog(Workspace=e2WS, LogName="vert_axis", LogType="String", LogText=label)
385
386	if unit[0] == 'Temperature':
387	nT = len(Tvalue)
388	if Tvalue[0] < Tvalue[nT-1]:
389	lo = 0
390	hi = nT-1
391	else:
392	lo = nT-1
393	hi = 0
394	text = 'Temperature range : '+str(Tvalue[lo])+' to '+str(Tvalue[hi])
395	AddSampleLog(Workspace=e1WS, LogName="temperature_normalise", LogType="String", LogText=str(Normalise))
396	if Normalise:
397	eltWS = ename+'_elt'
398	nhist = mtd[elfWS].getNumberHistograms()
399	for n in range(0,nhist):
400	ExtractSingleSpectrum(InputWorkspace=elfWS, OutputWorkspace='__lowTemp', WorkspaceIndex=n)
401	y = mtd['__lowTemp'].readY(0)
402	scale = 1.0/y[0]
403	Scale(InputWorkspace='__lowTemp', OutputWorkspace='__lowTemp', Factor=scale, Operation='Multiply')
404	if n == 0:
405	RenameWorkspace(InputWorkspace='__lowTemp', OutputWorkspace=eltWS)
406	else:
407	ConjoinWorkspaces(InputWorkspace1=eltWS, InputWorkspace2='__lowTemp')
408	# need to create a vertical axis here eg VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q1
409	CopyLogs(InputWorkspace='run1', OutputWorkspace=eltWS)
410	unitx = mtd[eltWS].getAxis(0).setUnit("Label")
411	unitx.setLabel(unit[0], unit[1])
412	addElwinLogs(elfWS, label, eRange, Range2)
413	DeleteWorkspace('run1')
414
415	unity = mtd[e1WS].getAxis(1).setUnit("Label")
416	unity.setLabel(unit[0], unit[1])
417	label = unit[0]+' / '+unit[1]
418	addElwinLogs(e1WS, label, eRange, Range2)
419
420	unity = mtd[e2WS].getAxis(1).setUnit("Label")
421	unity.setLabel(unit[0], unit[1])
422	addElwinLogs(e2WS, label, eRange, Range2)
423
424	unitx = mtd[elfWS].getAxis(0).setUnit("Label")
425	unitx.setLabel(unit[0], unit[1])
426	addElwinLogs(elfWS, label, eRange, Range2)
427
428	if Save:
429	e1_path = os.path.join(workdir, e1WS+'.nxs') # path name for nxs file
430	e2_path = os.path.join(workdir, e2WS+'.nxs') # path name for nxs file
431	elf_path = os.path.join(workdir, elfWS+'.nxs') # path name for nxs file
432
433	if Verbose:
434	logger.notice('Creating file : '+e1_path)
435	logger.notice('Creating file : '+e2_path)
436	logger.notice('Creating file : '+elf_path)
437
438	SaveNexusProcessed(InputWorkspace=e1WS, Filename=e1_path)
439	SaveNexusProcessed(InputWorkspace=e2WS, Filename=e2_path)
440	SaveNexusProcessed(InputWorkspace=elfWS, Filename=elf_path)
441	if Normalise:
442	elt_path = os.path.join(workdir, eltWS+'.nxs') # path name for nxs file
443	SaveNexusProcessed(InputWorkspace=eltWS, Filename=elt_path)
444	if Verbose:
445	logger.notice('Creating file : '+elt_path)
446
447	if Plot:
448	elwinPlot(label,e1WS,e2WS,elfWS)
449
450	EndTime('Elwin')
451	return e1WS,e2WS
452
453	# Add sample log to each of the workspaces created by Elwin
454	def addElwinLogs(ws, label, eRange, Range2):
455
456	AddSampleLog(Workspace=ws, LogName="vert_axis", LogType="String", LogText=label)
457	AddSampleLog(Workspace=ws, LogName="range1_start", LogType="Number", LogText=str(eRange[0]))
458	AddSampleLog(Workspace=ws, LogName="range1_end", LogType="Number", LogText=str(eRange[1]))
459	AddSampleLog(Workspace=ws, LogName="two_ranges", LogType="String", LogText=str(Range2))
460
461	if Range2:
462	AddSampleLog(Workspace=ws, LogName="range2_start", LogType="Number", LogText=str(eRange[2]))
463	AddSampleLog(Workspace=ws, LogName="range2_end", LogType="Number", LogText=str(eRange[3]))
464
465	#Plot each of the workspace output by elwin
466	def elwinPlot(label,eq1,eq2,elf):
467	plotElwinWorkspace(eq1, yAxisTitle='Elastic Intensity', setScale=True)
468	plotElwinWorkspace(eq2, yAxisTitle='log(Elastic Intensity)', setScale=True)
469	plotElwinWorkspace(elf, xAxisTitle=label)
470
471	#Plot a workspace generated by Elwin
472	def plotElwinWorkspace(ws, xAxisTitle=None, yAxisTitle=None, setScale=False):
473	ws = mtd[ws]
474	nBins = ws.blocksize()
475	lastX = ws.readX(0)[nBins-1]
476
477	nhist = ws.getNumberHistograms()
478
479	try:
480	graph = mp.plotSpectrum(ws, range(0,nhist))
481	except RuntimeError, e:
482	#User clicked cancel on plot so don't do anything
483	return None
484
485	layer = graph.activeLayer()
486
487	#set the x scale of the layer
488	if setScale:
489	layer.setScale(mp.Layer.Bottom, 0.0, lastX)
490
491	#set the title on the x-axis
492	if xAxisTitle:
493	layer.setAxisTitle(mp.Layer.Bottom, xAxisTitle)
494
495	#set the title on the y-axis
496	if yAxisTitle:
497	layer.setAxisTitle(mp.Layer.Left, yAxisTitle)
498
499	##############################################################################
500	# Fury
501	##############################################################################
502
503	def furyPlot(inWS, spec):
504	graph = mp.plotSpectrum(inWS, spec)
505	layer = graph.activeLayer()
506	layer.setScale(mp.Layer.Left, 0, 1.0)
507
508	def fury(samWorkspaces, res_file, rebinParam, RES=True, Save=False, Verbose=False,
509	Plot=False):
510	StartTime('Fury')
511	workdir = config['defaultsave.directory']
512	samTemp = samWorkspaces[0]
513	nsam,npt = CheckHistZero(samTemp)
514	Xin = mtd[samTemp].readX(0)
515	d1 = Xin[1]-Xin[0]
516	if d1 < 1e-8:
517	error = 'Data energy bin is zero'
518	logger.notice('ERROR *** ' + error)
519	sys.exit(error)
520	d2 = Xin[npt-1]-Xin[npt-2]
521	dmin = min(d1,d2)
522	pars = rebinParam.split(',')
523	if (float(pars[1]) <= dmin):
524	error = 'EWidth = ' + pars[1] + ' < smallest Eincr = ' + str(dmin)
525	logger.notice('ERROR *** ' + error)
526	sys.exit(error)
527	outWSlist = []
528	# Process RES Data Only Once
529	if Verbose:
530	logger.notice('Reading RES file : '+res_file)
531	LoadNexus(Filename=res_file, OutputWorkspace='res_data') # RES
532	CheckAnalysers(samTemp,'res_data',Verbose)
533	nres,nptr = CheckHistZero('res_data')
534	if nres > 1:
535	CheckHistSame(samTemp,'Sample','res_data','Resolution')
536	Rebin(InputWorkspace='res_data', OutputWorkspace='res_data', Params=rebinParam)
537	Integration(InputWorkspace='res_data', OutputWorkspace='res_int')
538	ConvertToPointData(InputWorkspace='res_data', OutputWorkspace='res_data')
539	ExtractFFTSpectrum(InputWorkspace='res_data', OutputWorkspace='res_fft', FFTPart=2)
540	Divide(LHSWorkspace='res_fft', RHSWorkspace='res_int', OutputWorkspace='res')
541	for samWs in samWorkspaces:
542	(direct, filename) = os.path.split(samWs)
543	(root, ext) = os.path.splitext(filename)
544	Rebin(InputWorkspace=samWs, OutputWorkspace='sam_data', Params=rebinParam)
545	Integration(InputWorkspace='sam_data', OutputWorkspace='sam_int')
546	ConvertToPointData(InputWorkspace='sam_data', OutputWorkspace='sam_data')
547	ExtractFFTSpectrum(InputWorkspace='sam_data', OutputWorkspace='sam_fft', FFTPart=2)
548	Divide(LHSWorkspace='sam_fft', RHSWorkspace='sam_int', OutputWorkspace='sam')
549	# Create save file name
550	savefile = getWSprefix('sam_data', root) + 'iqt'
551	outWSlist.append(savefile)
552	Divide(LHSWorkspace='sam', RHSWorkspace='res', OutputWorkspace=savefile)
553	#Cleanup Sample Files
554	DeleteWorkspace('sam_data')
555	DeleteWorkspace('sam_int')
556	DeleteWorkspace('sam_fft')
557	DeleteWorkspace('sam')
558	# Crop nonsense values off workspace
559	bin = int(math.ceil(mtd[savefile].blocksize()/2.0))
560	binV = mtd[savefile].dataX(0)[bin]
561	CropWorkspace(InputWorkspace=savefile, OutputWorkspace=savefile, XMax=binV)
562	if Save:
563	opath = os.path.join(workdir, savefile+'.nxs') # path name for nxs file
564	SaveNexusProcessed(InputWorkspace=savefile, Filename=opath)
565	if Verbose:
566	logger.notice('Output file : '+opath)
567	# Clean Up RES files
568	DeleteWorkspace('res_data')
569	DeleteWorkspace('res_int')
570	DeleteWorkspace('res_fft')
571	DeleteWorkspace('res')
572	if Plot:
573	specrange = range(0,mtd[outWSlist[0]].getNumberHistograms())
574	furyPlot(outWSlist, specrange)
575	EndTime('Fury')
576	return outWSlist
577
578	##############################################################################
579	# FuryFit
580	##############################################################################
581
582	def getFuryFitOption(option):
583	nopt = len(option)
584	if nopt == 2:
585	npeak = option[0]
586	type = option[1]
587	elif nopt == 4:
588	npeak = '2'
589	type = 'SE'
590	else:
591	error = 'Bad option : ' +option
592	logger.notice('ERROR *** ' + error)
593	sys.exit(error)
594	return npeak, type
595
596	def furyfitParsToWS(Table, Data, option):
597	npeak, type = getFuryFitOption(option)
598	Q = createQaxis(Data)
599	nQ = len(Q)
600	ws = mtd[Table]
601	rCount = ws.rowCount()
602	cCount = ws.columnCount()
603	cName = ws.getColumnNames()
604	Qa = np.array(Q)
605	A0v = ws.column(1) #bgd value
606	A0e = ws.column(2) #bgd error
607	Iy1 = ws.column(5) #intensity1 value
608	Ie1 = ws.column(2) #intensity1 error = bgd
609	dataX = Qa
610	dataY = np.array(A0v)
611	dataE = np.array(A0e)
612	names = cName[1]
613	dataX = np.append(dataX,Qa)
614	dataY = np.append(dataY,np.array(Iy1))
615	dataE = np.append(dataE,np.array(Ie1))
616	names += ","+cName[5]
617	Ty1 = ws.column(7) #tau1 value
618	Te1 = ws.column(8) #tau1 error
619	dataX = np.append(dataX,Qa)
620	dataY = np.append(dataY,np.array(Ty1))
621	dataE = np.append(dataE,np.array(Te1))
622	names += ","+cName[7]
623	nSpec = 3
624	if npeak == '1' and type == 'S':
625	By1 = ws.column(9) #beta1 value
626	Be1 = ws.column(10) #beta2 error
627	dataX = np.append(dataX,Qa)
628	dataY = np.append(dataY,np.array(By1))
629	dataE = np.append(dataE,np.array(Be1))
630	names += ","+cName[9]
631	nSpec += 1
632	if npeak == '2':
633	Iy2 = ws.column(9) #intensity2 value
634	Ie2 = ws.column(10) #intensity2 error
635	dataX = np.append(dataX,Qa)
636	dataY = np.append(dataY,np.array(Iy2))
637	dataE = np.append(dataE,np.array(Ie2))
638	names += ","+cName[9]
639	nSpec += 1
640	Ty2 = ws.column(11) #tau2 value
641	Te2 = ws.column(12) #tau2 error
642	dataX = np.append(dataX,Qa)
643	dataY = np.append(dataY,np.array(Ty2))
644	dataE = np.append(dataE,np.array(Te2))
645	names += ","+cName[11]
646	nSpec += 1
647	wsname = Table + "_Workspace"
648	CreateWorkspace(OutputWorkspace=wsname, DataX=dataX, DataY=dataY, DataE=dataE,
649	Nspec=nSpec, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
650	VerticalAxisValues=names)
651	return wsname
652
653	def createFurySeqResFun(ties, par, option):
654	npeak, type = getFuryFitOption(option)
655	fun = 'name=LinearBackground,A0='+str(par[0])+',A1=0,ties=(A1=0);'
656	if npeak == '1' and type == 'E':
657	fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)),Intensity='+str(par[1])+',Tau='+str(par[2])
658	if npeak == '1' and type == 'S':
659	fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity='+str(par[1])+',Tau='+str(par[2])+',Beta='+str(par[3])
660	if ties:
661	fun += ';ties=(f1.Intensity=1-f0.A0)'
662	return fun
663
664	def getFurySeqResult(inputWS, outNm, option, Verbose):
665	logger.notice('Option : ' +option)
666	npeak, type = getFuryFitOption(option)
667	params = mtd[outNm+'_Parameters']
668	A0 = params.column(1) #bgd value
669	I1 = params.column(5) #intensity1 value
670	T1 = params.column(7) #tau1 value
671	if npeak == '1' and type == 'S':
672	B1 = params.column(9) #beta1 value
673	if npeak == '2':
674	I2 = params.column(9) #intensity2 value
675	T2 = params.column(11) #tau2 value
676	nHist = mtd[inputWS].getNumberHistograms()
677	for i in range(nHist):
678	paras = [A0[i], I1[i], T1[i]]
679	if npeak == '1' and type == 'S':
680	paras.append(B1[i])
681	if npeak == '2':
682	paras.append(I2[i])
683	paras.append(T2[i])
684	func = createFurySeqResFun(True, paras, option)
685	if Verbose:
686	logger.notice('Fit func : '+func)
687	fitWS = outNm + '_Result_'
688	fout = fitWS + str(i)
689	Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i,Output=fout,MaxIterations=0)
690	unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
691	unitx.setLabel('Time' , 'ns')
692	RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
693	DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
694	DeleteWorkspace(fitWS+str(i)+'_Parameters')
695	if i == 0:
696	group = fout
697	else:
698	group += ',' + fout
699	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
700
701	def furyfitPlotSeq(inputWS, Plot):
702	nHist = mtd[inputWS].getNumberHistograms()
703	if ( Plot == 'All' ):
704	mp.plotSpectrum(inputWS, range(0, nHist), True)
705	return
706	plotSpecs = []
707	if ( Plot == 'Intensity' ):
708	res = 'Intensity$'
709	if ( Plot == 'Tau' ):
710	res = 'Tau$'
711	elif ( Plot == 'Beta' ):
712	res = 'Beta$'
713	for i in range(0, nHist):
714	title = mtd[inputWS].getAxis(1).label(i)
715	if ( re.search(res, title) ):
716	plotSpecs.append(i)
717	mp.plotSpectrum(inputWS, plotSpecs, True)
718
719	def furyfitSeq(inputWS, func, ftype, startx, endx, Save, Plot, Verbose=False):
720	StartTime('FuryFit')
721	workdir = config['defaultsave.directory']
722	input = inputWS+',i0'
723	nHist = mtd[inputWS].getNumberHistograms()
724	for i in range(1,nHist):
725	input += ';'+inputWS+',i'+str(i)
726	outNm = getWSprefix(inputWS) + 'fury_' + ftype + "0_to_" + str(nHist-1)
727	option = ftype[:-2]
728	if Verbose:
729	logger.notice('Option: '+option)
730	logger.notice(func)
731	PlotPeakByLogValue(Input=input, OutputWorkspace=outNm, Function=func,
732	StartX=startx, EndX=endx, FitType='Sequential')
733	fitWS = furyfitParsToWS(outNm, inputWS, option)
734	RenameWorkspace(InputWorkspace=outNm, OutputWorkspace=outNm+"_Parameters")
735	CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=inputWS, XMin=startx, XMax=endx)
736	getFurySeqResult(inputWS, outNm, option, Verbose)
737	if Save:
738	opath = os.path.join(workdir, fitWS+'.nxs') # path name for nxs file
739	SaveNexusProcessed(InputWorkspace=fitWS, Filename=opath)
740	if Verbose:
741	logger.notice('Output file : '+opath)
742	if ( Plot != 'None' ):
743	furyfitPlotSeq(fitWS, Plot)
744	EndTime('FuryFit')
745	return mtd[fitWS]
746
747	def furyfitMultParsToWS(Table, Data):
748	# Q = createQaxis(Data)
749	theta,Q = GetThetaQ(Data)
750	ws = mtd[Table+'_Parameters']
751	rCount = ws.rowCount()
752	cCount = ws.columnCount()
753	nSpec = ( rCount - 1 ) / 5
754	val = ws.column(1) #value
755	err = ws.column(2) #error
756	dataX = []
757	A0val = []
758	A0err = []
759	Ival = []
760	Ierr = []
761	Tval = []
762	Terr = []
763	Bval = []
764	Berr = []
765	for spec in range(0,nSpec):
766	n1 = spec*5
767	A0 = n1
768	A1 = n1+1
769	int = n1+2 #intensity value
770	tau = n1+3 #tau value
771	beta = n1+4 #beta value
772	dataX.append(Q[spec])
773	A0val.append(val[A0])
774	A0err.append(err[A0])
775	Ival.append(val[int])
776	Ierr.append(err[int])
777	Tval.append(val[tau])
778	Terr.append(err[tau])
779	Bval.append(val[beta])
780	Berr.append(err[beta])
781	nQ = len(dataX)
782	Qa = np.array(dataX)
783	dataY = np.array(A0val)
784	dataE = np.array(A0err)
785	dataY = np.append(dataY,np.array(Ival))
786	dataE = np.append(dataE,np.array(Ierr))
787	dataY = np.append(dataY,np.array(Tval))
788	dataE = np.append(dataE,np.array(Terr))
789	dataY = np.append(dataY,np.array(Bval))
790	dataE = np.append(dataE,np.array(Berr))
791	names = 'A0,Intensity,Tau,Beta'
792	suffix = 'Workspace'
793	wsname = Table + '_' + suffix
794	CreateWorkspace(OutputWorkspace=wsname, DataX=Qa, DataY=dataY, DataE=dataE,
795	Nspec=4, UnitX='MomentumTransfer', VerticalAxisUnit='Text',
796	VerticalAxisValues=names)
797	return wsname
798
799	def furyfitPlotMult(inputWS, Plot):
800	nHist = mtd[inputWS].getNumberHistograms()
801	if ( Plot == 'All' ):
802	mp.plotSpectrum(inputWS, range(0, nHist))
803	return
804	plotSpecs = []
805	if ( Plot == 'Intensity' ):
806	mp.plotSpectrum(inputWS, 1, True)
807	if ( Plot == 'Tau' ):
808	mp.plotSpectrum(inputWS, 2, True)
809	elif ( Plot == 'Beta' ):
810	mp.plotSpectrum(inputWS, 3, True)
811
812
813	def createFuryMultFun(ties = True, function = ''):
814	fun = '(composite=CompositeFunction,$domains=i;'
815	fun += function
816	if ties:
817	fun += ';ties=(f1.Intensity=1-f0.A0)'
818	fun += ');'
819	return fun
820
821	def createFuryMultResFun(ties = True, A0 = 0.02, Intensity = 0.98 ,Tau = 0.025, Beta = 0.8):
822	fun = '(composite=CompositeFunction,$domains=i;'
823	fun += 'name=LinearBackground,A0='+str(A0)+',A1=0,ties=(A1=0);'
824	fun += 'name=UserFunction,Formula=Intensity*exp(-(x/Tau)^Beta),Intensity='+str(Intensity)+',Tau='+str(Tau)+',Beta='+str(Beta)
825	if ties:
826	fun += ';ties=(f1.Intensity=1-f0.A0)'
827	fun += ');'
828	return fun
829
830	def getFuryMultResult(inputWS, outNm, function, Verbose):
831	params = mtd[outNm+'_Parameters']
832	nHist = mtd[inputWS].getNumberHistograms()
833	for i in range(nHist):
834	j = 5 * i
835	# assert( params.row(j)['Name'][3:] == 'f0.A0' )
836	A0 = params.row(j)['Value']
837	A1 = params.row(j + 1)['Value']
838	Intensity = params.row(j + 2)['Value']
839	Tau = params.row(j + 3)['Value']
840	Beta = params.row(j + 4)['Value']
841	func = createFuryMultResFun(True, A0, Intensity ,Tau, Beta)
842	if Verbose:
843	logger.notice('Fit func : '+func)
844	fitWS = outNm + '_Result_'
845	fout = fitWS + str(i)
846	Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=i,Output=fout,MaxIterations=0)
847	unitx = mtd[fout+'_Workspace'].getAxis(0).setUnit("Label")
848	unitx.setLabel('Time' , 'ns')
849	RenameWorkspace(InputWorkspace=fout+'_Workspace', OutputWorkspace=fout)
850	DeleteWorkspace(fitWS+str(i)+'_NormalisedCovarianceMatrix')
851	DeleteWorkspace(fitWS+str(i)+'_Parameters')
852	if i == 0:
853	group = fout
854	else:
855	group += ',' + fout
856	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fitWS[:-1])
857
858	def furyfitMult(inputWS, function, ftype, startx, endx, Save, Plot, Verbose=False):
859	StartTime('FuryFit Mult')
860	workdir = config['defaultsave.directory']
861	option = ftype[:-2]
862	if Verbose:
863	logger.notice('Option: '+option)
864	logger.notice('Function: '+function)
865	nHist = mtd[inputWS].getNumberHistograms()
866	outNm = inputWS[:-3] + 'fury_mult'
867	f1 = createFuryMultFun(True, function)
868	func= 'composite=MultiDomainFunction,NumDeriv=1;'
869	ties='ties=('
870	kwargs = {}
871	for i in range(0,nHist):
872	func+=f1
873	if i > 0:
874	ties += 'f' + str(i) + '.f1.Beta=f0.f1.Beta'
875	if i < nHist-1:
876	ties += ','
877	kwargs['InputWorkspace_' + str(i)] = inputWS
878	kwargs['WorkspaceIndex_' + str(i)] = i
879	ties+=')'
880	func += ties
881	CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=inputWS, XMin=startx, XMax=endx)
882	Fit(Function=func,InputWorkspace=inputWS,WorkspaceIndex=0,Output=outNm,**kwargs)
883	outWS = furyfitMultParsToWS(outNm, inputWS)
884	getFuryMultResult(inputWS, outNm, function, Verbose)
885	if Save:
886	opath = os.path.join(workdir, outWS+'.nxs') # path name for nxs file
887	SaveNexusProcessed(InputWorkspace=outWS, Filename=opath)
888	rpath = os.path.join(workdir, outNm+'_result.nxs') # path name for nxs file
889	SaveNexusProcessed(InputWorkspace=outNm+'_result', Filename=rpath)
890	if Verbose:
891	logger.notice('Output file : '+opath)
892	logger.notice('Output file : '+rpath)
893	if ( Plot != 'None' ):
894	furyfitPlotMult(outWS, Plot)
895	EndTime('FuryFit')
896
897	##############################################################################
898	# MSDFit
899	##############################################################################
900
901	def msdfitParsToWS(Table, xData):
902	dataX = xData
903	ws = mtd[Table+'_Table']
904	rCount = ws.rowCount()
905	yA0 = ws.column(1)
906	eA0 = ws.column(2)
907	yA1 = ws.column(3)
908	dataY1 = map(lambda x : -x, yA1)
909	eA1 = ws.column(4)
910	wsname = Table
911	CreateWorkspace(OutputWorkspace=wsname+'_a0', DataX=dataX, DataY=yA0, DataE=eA0,
912	Nspec=1, UnitX='')
913	CreateWorkspace(OutputWorkspace=wsname+'_a1', DataX=dataX, DataY=dataY1, DataE=eA1,
914	Nspec=1, UnitX='')
915	group = wsname+'_a0,'+wsname+'_a1'
916	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=wsname)
917	return wsname
918
919	def msdfitPlotSeq(inputWS, xlabel):
920	msd_plot = mp.plotSpectrum(inputWS+'_a1',0,True)
921	msd_layer = msd_plot.activeLayer()
922	msd_layer.setAxisTitle(mp.Layer.Bottom,xlabel)
923	msd_layer.setAxisTitle(mp.Layer.Left,'<u2>')
924
925	def msdfitPlotFits(calcWS, n):
926	mfit_plot = mp.plotSpectrum(calcWS+'_0',[0,1],True)
927	mfit_layer = mfit_plot.activeLayer()
928	mfit_layer.setAxisTitle(mp.Layer.Left,'log(Elastic Intensity)')
929
930	def msdfit(inputs, startX, endX, Save=False, Verbose=False, Plot=True):
931	StartTime('msdFit')
932	workdir = config['defaultsave.directory']
933	log_type = 'sample'
934	file = inputs[0]
935	(direct, filename) = os.path.split(file)
936	(root, ext) = os.path.splitext(filename)
937	(instr, first) = getInstrRun(filename)
938	if Verbose:
939	logger.notice('Reading Run : '+file)
940	LoadNexusProcessed(FileName=file, OutputWorkspace=root)
941	nHist = mtd[root].getNumberHistograms()
942	file_list = []
943	run_list = []
944	ws = mtd[root]
945	ws_run = ws.getRun()
946	vertAxisValues = ws.getAxis(1).extractValues()
947	x_list = vertAxisValues
948	if 'vert_axis' in ws_run:
949	xlabel = ws_run.getLogData('vert_axis').value
950	for nr in range(0, nHist):
951	nsam,ntc = CheckHistZero(root)
952	lnWS = '__lnI_'+str(nr)
953	file_list.append(lnWS)
954	ExtractSingleSpectrum(InputWorkspace=root, OutputWorkspace=lnWS,
955	WorkspaceIndex=nr)
956	if (nr == 0):
957	run_list = lnWS
958	else:
959	run_list += ';'+lnWS
960	mname = root[:-4]
961	msdWS = mname+'_msd'
962	if Verbose:
963	logger.notice('Fitting Runs '+mname)
964	logger.notice('Q-range from '+str(startX)+' to '+str(endX))
965	function = 'name=LinearBackground, A0=0, A1=0'
966	PlotPeakByLogValue(Input=run_list, OutputWorkspace=msdWS+'_Table', Function=function,
967	StartX=startX, EndX=endX, FitType = 'Sequential')
968	msdfitParsToWS(msdWS, x_list)
969	nr = 0
970	fitWS = mname+'_Fit'
971	calcWS = mname+'_msd_Result'
972	a0 = mtd[msdWS+'_a0'].readY(0)
973	a1 = mtd[msdWS+'_a1'].readY(0)
974	for nr in range(0, nHist):
975	inWS = file_list[nr]
976	CropWorkspace(InputWorkspace=inWS,OutputWorkspace='__data',XMin=0.95startX,XMax=1.05endX)
977	dataX = mtd['__data'].readX(0)
978	nxd = len(dataX)
979	dataX = np.append(dataX,2*dataX[nxd-1]-dataX[nxd-2])
980	dataY = np.array(mtd['__data'].readY(0))
981	dataE = np.array(mtd['__data'].readE(0))
982	xd = []
983	yd = []
984	ed = []
985	for n in range(0,nxd):
986	line = a0[nr] - a1[nr]*dataX[n]
987	xd.append(dataX[n])
988	yd.append(line)
989	ed.append(0.0)
990	xd.append(dataX[nxd])
991	dataX = np.append(dataX,np.array(xd))
992	dataY = np.append(dataY,np.array(yd))
993	dataE = np.append(dataE,np.array(ed))
994	fout = calcWS +'_'+ str(nr)
995	CreateWorkspace(OutputWorkspace=fout, DataX=dataX, DataY=dataY, DataE=dataE,
996	Nspec=2, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues='Data,Calc')
997	if nr == 0:
998	gro = fout
999	else:
1000	gro += ',' + fout
1001	DeleteWorkspace(inWS)
1002	DeleteWorkspace('__data')
1003	GroupWorkspaces(InputWorkspaces=gro,OutputWorkspace=calcWS)
1004	CopyLogs(InputWorkspace=root, OutputWorkspace=msdWS)
1005	AddSampleLog(Workspace=msdWS, LogName="start_x", LogType="String", LogText=str(startX))
1006	AddSampleLog(Workspace=msdWS, LogName="end_x", LogType="String", LogText=str(endX))
1007	if Plot:
1008	msdfitPlotSeq(msdWS, xlabel)
1009	msdfitPlotFits(calcWS, 0)
1010	if Save:
1011	msd_path = os.path.join(workdir, msdWS+'.nxs') # path name for nxs file
1012	SaveNexusProcessed(InputWorkspace=msdWS, Filename=msd_path, Title=msdWS)
1013	if Verbose:
1014	logger.notice('Output msd file : '+msd_path)
1015	EndTime('msdFit')
1016	return msdWS
1017
1018	def plotInput(inputfiles,spectra=[]):
1019	OneSpectra = False
1020	if len(spectra) != 2:
1021	spectra = [spectra[0], spectra[0]]
1022	OneSpectra = True
1023	workspaces = []
1024	for file in inputfiles:
1025	root = LoadNexus(Filename=file)
1026	if not OneSpectra:
1027	GroupDetectors(root, root,
1028	DetectorList=range(spectra[0],spectra[1]+1) )
1029	workspaces.append(root)
1030	if len(workspaces) > 0:
1031	graph = mp.plotSpectrum(workspaces,0)
1032	layer = graph.activeLayer().setTitle(", ".join(workspaces))
1033
1034	##############################################################################
1035	# Corrections
1036	##############################################################################
1037
1038	def CubicFit(inputWS, spec, Verbose=False):
1039	'''Uses the Mantid Fit Algorithm to fit a quadratic to the inputWS
1040	parameter. Returns a list containing the fitted parameter values.'''
1041	function = 'name=Quadratic, A0=1, A1=0, A2=0'
1042	fit = Fit(Function=function, InputWorkspace=inputWS, WorkspaceIndex=spec,
1043	CreateOutput=True, Output='Fit')
1044	table = mtd['Fit_Parameters']
1045	A0 = table.cell(0,1)
1046	A1 = table.cell(1,1)
1047	A2 = table.cell(2,1)
1048	Abs = [A0, A1, A2]
1049	if Verbose:
1050	logger.notice('Group '+str(spec)+' of '+inputWS+' ; fit coefficients are : '+str(Abs))
1051	return Abs
1052
1053	def applyCorrections(inputWS, canWS, corr, Verbose=False):
1054	'''Through the PolynomialCorrection algorithm, makes corrections to the
1055	input workspace based on the supplied correction values.'''
1056	# Corrections are applied in Lambda (Wavelength)
1057	efixed = getEfixed(inputWS) # Get efixed
1058	theta,Q = GetThetaQ(inputWS)
1059	sam_name = getWSprefix(inputWS)
1060	ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='Wavelength',
1061	EMode='Indirect', EFixed=efixed)
1062
1063	nameStem = corr[:-4]
1064	if canWS != '':
1065	(instr, can_run) = getInstrRun(canWS)
1066	corrections = [nameStem+'_ass', nameStem+'_assc', nameStem+'_acsc', nameStem+'_acc']
1067	CorrectedWS = sam_name +'Correct_'+ can_run
1068	ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='Wavelength',
1069	EMode='Indirect', EFixed=efixed)
1070	else:
1071	corrections = [nameStem+'_ass']
1072	CorrectedWS = sam_name +'Corrected'
1073	nHist = mtd[inputWS].getNumberHistograms()
1074	# Check that number of histograms in each corrections workspace matches
1075	# that of the input (sample) workspace
1076	for ws in corrections:
1077	if ( mtd[ws].getNumberHistograms() != nHist ):
1078	raise ValueError('Mismatch: num of spectra in '+ws+' and inputWS')
1079	# Workspaces that hold intermediate results
1080	CorrectedSampleWS = '__csam'
1081	CorrectedCanWS = '__ccan'
1082	for i in range(0, nHist): # Loop through each spectra in the inputWS
1083	ExtractSingleSpectrum(InputWorkspace=inputWS, OutputWorkspace=CorrectedSampleWS,
1084	WorkspaceIndex=i)
1085	if ( len(corrections) == 1 ):
1086	Ass = CubicFit(corrections[0], i, Verbose)
1087	PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
1088	Coefficients=Ass, Operation='Divide')
1089	if ( i == 0 ):
1090	CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
1091	else:
1092	ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS)
1093	else:
1094	ExtractSingleSpectrum(InputWorkspace=canWS, OutputWorkspace=CorrectedCanWS,
1095	WorkspaceIndex=i)
1096	Acc = CubicFit(corrections[3], i, Verbose)
1097	PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
1098	Coefficients=Acc, Operation='Divide')
1099	Acsc = CubicFit(corrections[2], i, Verbose)
1100	PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS,
1101	Coefficients=Acsc, Operation='Multiply')
1102	Minus(LHSWorkspace=CorrectedSampleWS, RHSWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedSampleWS)
1103	Assc = CubicFit(corrections[1], i, Verbose)
1104	PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS,
1105	Coefficients=Assc, Operation='Divide')
1106	if ( i == 0 ):
1107	CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS)
1108	else:
1109	ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS,
1110	CheckOverlapping=False)
1111	ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='DeltaE',
1112	EMode='Indirect', EFixed=efixed)
1113	ConvertUnits(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS, Target='DeltaE',
1114	EMode='Indirect', EFixed=efixed)
1115	ConvertSpectrumAxis(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_rqw',
1116	Target='ElasticQ', EMode='Indirect', EFixed=efixed)
1117	RenameWorkspace(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_red')
1118	if canWS != '':
1119	DeleteWorkspace(CorrectedCanWS)
1120	ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='DeltaE',
1121	EMode='Indirect', EFixed=efixed)
1122	DeleteWorkspace('Fit_NormalisedCovarianceMatrix')
1123	DeleteWorkspace('Fit_Parameters')
1124	DeleteWorkspace('Fit_Workspace')
1125	DeleteWorkspace(corr)
1126	return CorrectedWS
1127
1128	def abscorFeeder(sample, container, geom, useCor, corrections, Verbose=False, ScaleOrNotToScale=False, factor=1, Save=False,
1129	PlotResult='None', PlotContrib=False):
1130	'''Load up the necessary files and then passes them into the main
1131	applyCorrections routine.'''
1132	StartTime('ApplyCorrections')
1133	workdir = config['defaultsave.directory']
1134	s_hist,sxlen = CheckHistZero(sample)
1135	sam_name = getWSprefix(sample)
1136	efixed = getEfixed(sample)
1137	if container != '':
1138	CheckAnalysers(sample,container,Verbose)
1139	CheckHistSame(sample,'Sample',container,'Container')
1140	(instr, can_run) = getInstrRun(container)
1141	if ScaleOrNotToScale:
1142	Scale(InputWorkspace=container, OutputWorkspace=container, Factor=factor, Operation='Multiply')
1143	if Verbose:
1144	logger.notice('Container scaled by '+str(factor))
1145	if useCor:
1146	if Verbose:
1147	text = 'Correcting sample ' + sample
1148	if container != '':
1149	text += ' with ' + container
1150	logger.notice(text)
1151
1152	cor_result = applyCorrections(sample, container, corrections, Verbose)
1153	rws = mtd[cor_result+'_red']
1154	outNm= cor_result + '_Result_'
1155
1156	if Save:
1157	cred_path = os.path.join(workdir,cor_result+'_red.nxs')
1158	SaveNexusProcessed(InputWorkspace=cor_result+'_red',Filename=cred_path)
1159	if Verbose:
1160	logger.notice('Output file created : '+cred_path)
1161	calc_plot = [cor_result+'_red',sample]
1162	res_plot = cor_result+'_rqw'
1163	else:
1164	if ( container == '' ):
1165	sys.exit('ERROR *** Invalid options - nothing to do!')
1166	else:
1167	sub_result = sam_name +'Subtract_'+ can_run
1168	if Verbose:
1169	logger.notice('Subtracting '+container+' from '+sample)
1170	Minus(LHSWorkspace=sample,RHSWorkspace=container,OutputWorkspace=sub_result)
1171	ConvertSpectrumAxis(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_rqw',
1172	Target='ElasticQ', EMode='Indirect', EFixed=efixed)
1173	RenameWorkspace(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_red')
1174	rws = mtd[sub_result+'_red']
1175	outNm= sub_result + '_Result_'
1176	if Save:
1177	sred_path = os.path.join(workdir,sub_result+'_red.nxs')
1178	SaveNexusProcessed(InputWorkspace=sub_result+'_red',Filename=sred_path)
1179	if Verbose:
1180	logger.notice('Output file created : '+sred_path)
1181	res_plot = sub_result+'_rqw'
1182	if (PlotResult != 'None'):
1183	plotCorrResult(res_plot,PlotResult)
1184	if ( container != '' ):
1185	sws = mtd[sample]
1186	cws = mtd[container]
1187	names = 'Sample,Can,Calc'
1188	for i in range(0, s_hist): # Loop through each spectra in the inputWS
1189	dataX = np.array(sws.readX(i))
1190	dataY = np.array(sws.readY(i))
1191	dataE = np.array(sws.readE(i))
1192	dataX = np.append(dataX,np.array(cws.readX(i)))
1193	dataY = np.append(dataY,np.array(cws.readY(i)))
1194	dataE = np.append(dataE,np.array(cws.readE(i)))
1195	dataX = np.append(dataX,np.array(rws.readX(i)))
1196	dataY = np.append(dataY,np.array(rws.readY(i)))
1197	dataE = np.append(dataE,np.array(rws.readE(i)))
1198	fout = outNm + str(i)
1199	CreateWorkspace(OutputWorkspace=fout, DataX=dataX, DataY=dataY, DataE=dataE,
1200	Nspec=3, UnitX='DeltaE', VerticalAxisUnit='Text', VerticalAxisValues=names)
1201	if i == 0:
1202	group = fout
1203	else:
1204	group += ',' + fout
1205	GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=outNm[:-1])
1206	if PlotContrib:
1207	plotCorrContrib(outNm+'0',[0,1,2])
1208	if Save:
1209	res_path = os.path.join(workdir,outNm[:-1]+'.nxs')
1210	SaveNexusProcessed(InputWorkspace=outNm[:-1],Filename=res_path)
1211	if Verbose:
1212	logger.notice('Output file created : '+res_path)
1213	EndTime('ApplyCorrections')
1214
1215	def plotCorrResult(inWS,PlotResult):
1216	nHist = mtd[inWS].getNumberHistograms()
1217	if (PlotResult == 'Spectrum' or PlotResult == 'Both'):
1218	if nHist >= 10: #only plot up to 10 hists
1219	nHist = 10
1220	plot_list = []
1221	for i in range(0, nHist):
1222	plot_list.append(i)
1223	res_plot=mp.plotSpectrum(inWS,plot_list)
1224	if (PlotResult == 'Contour' or PlotResult == 'Both'):
1225	if nHist >= 5: #needs at least 5 hists for a contour
1226	mp.importMatrixWorkspace(inWS).plotGraph2D()
1227
1228	def plotCorrContrib(plot_list,n):
1229	con_plot=mp.plotSpectrum(plot_list,n)

Download in other formats:

Original Format