Context Navigation

Back to Ticket #7835

Ticket #7835: IndirectDataAnalysis.py

File IndirectDataAnalysis.py, 47.0 KB (added by Samuel Jackson, 7 years ago)
Spencer's updated file with mods on lines 1091-2, at 1053 & around 1031-32

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

Download in other formats:

Original Format