Context Navigation

Back to Ticket #7977

Ticket #7977: IndirectDataAnalysis.py

File IndirectDataAnalysis.py, 48.8 KB (added by Samuel Jackson, 7 years ago)
IDA python script with updated Elwin functions

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

Download in other formats:

Original Format