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