From d190bd4056ea0b018bc7d88436aa730fdaad769c Mon Sep 17 00:00:00 2001 From: Dan Nixon Date: Fri, 15 Aug 2014 15:32:18 +0100 Subject: [PATCH 1/3] Revert "IndirectDataAnalysis refactoring" This reverts commit d498250fca8c42284530060d6b3f27ec914559a6. --- .../scripts/Inelastic/IndirectDataAnalysis.py | 886 +++++++++--------- 1 file changed, 445 insertions(+), 441 deletions(-) diff --git a/Code/Mantid/scripts/Inelastic/IndirectDataAnalysis.py b/Code/Mantid/scripts/Inelastic/IndirectDataAnalysis.py index 90ec0d4950c1..d405ab8c3306 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectDataAnalysis.py +++ b/Code/Mantid/scripts/Inelastic/IndirectDataAnalysis.py @@ -1,8 +1,8 @@ from IndirectImport import import_mantidplot -MTD_PLOT = import_mantidplot() +mp = import_mantidplot() from IndirectCommon import * -import math, re, sys, os.path, numpy as np +import math, re, os.path, numpy as np from mantid.simpleapi import * from mantid.api import TextAxis from mantid import * @@ -22,13 +22,13 @@ def segment(l, fromIndex, toIndex): def trimData(nSpec, vals, min, max): result = [] - chunk_size = len(vals) / nSpec + chunkSize = len(vals) / nSpec assert min >= 0, 'trimData: min is less then zero' - assert max <= chunk_size - 1, 'trimData: max is greater than the number of spectra' + assert max <= chunkSize - 1, 'trimData: max is greater than the number of spectra' assert min <= max, 'trimData: min is greater than max' - chunks = split(vals, chunk_size) + chunks = split(vals,chunkSize) for chunk in chunks: - seg = segment(chunk, min, max) + seg = segment(chunk,min,max) for val in seg: result.append(val) return result @@ -56,10 +56,10 @@ def calculateEISF(params_table): #calculate EISF and EISF error total = height_y+amp_y - eisf_y = height_y/total + EISF_y = height_y/total total_error = height_e**2 + np.asarray(amp_e)**2 - eisf_e = eisf_y * np.sqrt((height_e**2/height_y**2) + (total_error/total**2)) + EISF_e = EISF_y * np.sqrt((height_e**2/height_y**2) + (total_error/total**2)) #append the calculated values to the table workspace col_name = amp_name[:-len('Amplitude')] + 'EISF' @@ -68,7 +68,7 @@ def calculateEISF(params_table): mtd[params_table].addColumn('double', col_name) mtd[params_table].addColumn('double', error_col_name) - for i, (value, error) in enumerate(zip(eisf_y, eisf_e)): + for i, (value, error) in enumerate(zip(EISF_y, EISF_e)): mtd[params_table].setCell(col_name, i, value) mtd[params_table].setCell(error_col_name, i, error) @@ -136,8 +136,8 @@ def confitSeq(inputWS, func, startX, endX, ftype, bgd, temperature=None, specMin RenameWorkspace(InputWorkspace=output_workspace, OutputWorkspace=output_workspace + "_Parameters") fit_workspaces = mtd[output_workspace + '_Workspaces'].getNames() - for i, workspace in enumerate(fit_workspaces): - RenameWorkspace(workspace, OutputWorkspace=output_workspace + '_' + str(i+specMin) + '_Workspace') + for i, ws in enumerate(fit_workspaces): + RenameWorkspace(ws, OutputWorkspace=output_workspace + '_' + str(i+specMin) + '_Workspace') if Save: # path name for nxs file @@ -164,7 +164,7 @@ def GetTemperature(root, tempWS, log_type, Verbose): facility = config.getFacility() pad_num = facility.instrument(instr).zeroPadding(int(run_number)) zero_padding = '0' * (pad_num - len(run_number)) - + run_name = instr + zero_padding + run_number log_name = run_name.upper() + '.log' @@ -174,16 +174,16 @@ def GetTemperature(root, tempWS, log_type, Verbose): # test logs in WS tmp = run[log_type].value temp = tmp[len(tmp)-1] - + if Verbose: mess = ' Run : '+run_name +' ; Temperature in log = '+str(temp) logger.notice(mess) - else: + else: # logs not in WS logger.warning('Log parameter not found in workspace. Searching for log file.') log_path = FileFinder.getFullPath(log_name) - - if log_path != '': + + if log_path != '': # get temperature from log file LoadLog(Workspace=tempWS, Filename=log_path) run_logs = mtd[tempWS].getRun() @@ -192,46 +192,46 @@ def GetTemperature(root, tempWS, log_type, Verbose): mess = ' Run : '+run_name+' ; Temperature in file = '+str(temp) logger.warning(mess) else: - # can't find log file + # can't find log file temp = int(run_name[-3:]) unit = ['Run-number', 'last 3 digits'] mess = ' Run : '+run_name +' ; Temperature file not found' logger.warning(mess) - return temp, unit + return temp,unit -def elwin(inputFiles, eRange, log_type='sample', Normalise=False, - Save=False, Verbose=False, Plot=False): +def elwin(inputFiles, eRange, log_type='sample', Normalise = False, + Save=False, Verbose=False, Plot=False): StartTime('ElWin') workdir = config['defaultsave.directory'] - CheckXrange(eRange, 'Energy') - temp_ws = '__temp' - range2_a = (len(eRange) == 4) + CheckXrange(eRange,'Energy') + tempWS = '__temp' + Range2 = ( len(eRange) == 4 ) if Verbose: range1 = str(eRange[0])+' to '+str(eRange[1]) - if len(eRange) == 4: + if ( len(eRange) == 4 ): range2 = str(eRange[2])+' to '+str(eRange[3]) logger.notice('Using 2 energy ranges from '+range1+' & '+range2) - elif len(eRange) == 2: + elif ( len(eRange) == 2 ): logger.notice('Using 1 energy range from '+range1) nr = 0 - input_runs = sorted(inputFiles) - for run_file in input_runs: - (_, file_name) = os.path.split(run_file) - (root, _) = os.path.splitext(file_name) - LoadNexus(Filename=file, OutputWorkspace=temp_ws) - CheckHistZero(temp_ws) - (xval, unit) = GetTemperature(root, temp_ws, log_type, Verbose) + inputRuns = sorted(inputFiles) + for file in inputRuns: + (direct, file_name) = os.path.split(file) + (root, ext) = os.path.splitext(file_name) + LoadNexus(Filename=file, OutputWorkspace=tempWS) + nsam,ntc = CheckHistZero(tempWS) + (xval, unit) = GetTemperature(root,tempWS,log_type, Verbose) if Verbose: - logger.notice('Reading file : '+run_file) - if len(eRange) == 4: - ElasticWindow(InputWorkspace=temp_ws, Range1Start=eRange[0], Range1End=eRange[1], + logger.notice('Reading file : '+file) + if ( len(eRange) == 4 ): + ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1], Range2Start=eRange[2], Range2End=eRange[3], OutputInQ='__eq1', OutputInQSquared='__eq2') - elif len(eRange) == 2: - ElasticWindow(InputWorkspace=temp_ws, Range1Start=eRange[0], Range1End=eRange[1], + elif ( len(eRange) == 2 ): + ElasticWindow(InputWorkspace=tempWS, Range1Start=eRange[0], Range1End=eRange[1], OutputInQ='__eq1', OutputInQSquared='__eq2') - (_, last) = getInstrRun(temp_ws) + (instr, last) = getInstrRun(tempWS) q1 = np.array(mtd['__eq1'].readX(0)) i1 = np.array(mtd['__eq1'].readY(0)) e1 = np.array(mtd['__eq1'].readE(0)) @@ -239,160 +239,161 @@ def elwin(inputFiles, eRange, log_type='sample', Normalise=False, q2 = np.array(mtd['__eq2'].readX(0)) i2 = np.array(mtd['__eq2'].readY(0)) e2 = np.array(mtd['__eq2'].readE(0)) - if nr == 0: + if (nr == 0): CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elf') - first = getWSprefix(temp_ws) - data_x1 = q1 - data_y1 = i1 - data_e1 = e1 - data_x2 = q2 - data_y2 = i2 - data_e2 = e2 - t_value = [xval] - t_error = [0.0] - t_axis = str(xval) + first = getWSprefix(tempWS) + datX1 = q1 + datY1 = i1 + datE1 = e1 + datX2 = q2 + datY2 = i2 + datE2 = e2 + Tvalue = [xval] + Terror = [0.0] + Taxis = str(xval) else: CloneWorkspace(InputWorkspace='__eq1', OutputWorkspace='__elftmp') ConjoinWorkspaces(InputWorkspace1='__elf', InputWorkspace2='__elftmp', CheckOverlapping=False) - data_x1 = np.append(data_x1, q1) - data_y1 = np.append(data_y1, i1) - data_e1 = np.append(data_e1, e1) - data_x2 = np.append(data_x2, q2) - data_y2 = np.append(data_y2, i2) - data_e2 = np.append(data_e2, e2) - t_value.append(xval) - t_error.append(0.0) - t_axis += ','+str(xval) + datX1 = np.append(datX1,q1) + datY1 = np.append(datY1,i1) + datE1 = np.append(datE1,e1) + datX2 = np.append(datX2,q2) + datY2 = np.append(datY2,i2) + datE2 = np.append(datE2,e2) + Tvalue.append(xval) + Terror.append(0.0) + Taxis += ','+str(xval) nr += 1 - txa = np.array(t_value) - num_q = len(q1) - for nq in range(0, num_q): + Txa = np.array(Tvalue) + Tea = np.array(Terror) + nQ = len(q1) + for nq in range(0,nQ): iq = [] eq = [] - for nt in range(0, len(t_value)): + for nt in range(0,len(Tvalue)): ii = mtd['__elf'].readY(nt) iq.append(ii[nq]) ie = mtd['__elf'].readE(nt) eq.append(ie[nq]) iqa = np.array(iq) eqa = np.array(eq) - if nq == 0: - data_tx = txa - data_ty = iqa - data_te = eqa + if (nq == 0): + datTx = Txa + datTy = iqa + datTe = eqa else: - data_tx = np.append(data_tx, txa) - data_ty = np.append(data_ty, iqa) - data_te = np.append(data_te, eqa) + datTx = np.append(datTx,Txa) + datTy = np.append(datTy,iqa) + datTe = np.append(datTe,eqa) DeleteWorkspace('__eq1') DeleteWorkspace('__eq2') DeleteWorkspace('__elf') - if nr == 1: + if (nr == 1): ename = first[:-1] else: ename = first+'to_'+last #check if temp was increasing or decreasing - if data_tx[0] > data_tx[-1]: + if(datTx[0] > datTx[-1]): # if so reverse data to follow natural ordering - data_tx = data_tx[::-1] - data_ty = data_ty[::-1] - data_te = data_te[::-1] + datTx = datTx[::-1] + datTy = datTy[::-1] + datTe = datTe[::-1] - elf_ws = ename + '_elf' - e1_ws = ename + '_eq1' - e2_ws = ename + '_eq2' + elfWS = ename+'_elf' + e1WS = ename+'_eq1' + e2WS = ename+'_eq2' #elt only created if we normalise - elt_ws = None - - wsnames = [elf_ws, e1_ws, e2_ws] + eltWS = None + wsnames = [elfWS, e1WS, e2WS] + #x,y,e data for the elf, e1 and e2 workspaces - data = [[data_tx, data_ty, data_te], - [data_x1, data_y1, data_e1], - [data_x2, data_y2, data_e2]] - + data = [[datTx, datTy, datTe], + [datX1, datY1, datE1], + [datX2, datY2, datE2]] + #x and vertical units for the elf, e1 and e2 workspaces xunits = ['Energy', 'MomentumTransfer', 'QSquared'] vunits = ['MomentumTransfer', 'Energy', 'Energy'] #vertical axis values for the elf, e1 and e2 workspaces - vvalues = [q1, t_axis, t_axis] - + vvalues = [q1, Taxis, Taxis] + #number of spectra in each workspace - nspecs = [num_q, nr, nr] + nspecs = [nQ, nr, nr] #x-axis units label label = unit[0]+' / '+unit[1] - ws_info = zip(wsnames, data, xunits, vunits, vvalues, nspecs) + wsInfo = zip(wsnames,data, xunits, vunits, vvalues, nspecs) #Create output workspaces and add sample logs - for wsname, wsdata, xunit, vunit, vvalue, nspec in ws_info: - x_val, y_val, e_val = wsdata + for wsname, wsdata, xunit, vunit, vvalue, nspec in wsInfo: + x, y, e = wsdata - CreateWorkspace(OutputWorkspace=wsname, DataX=x_val, DataY=y_val, DataE=e_val, + CreateWorkspace(OutputWorkspace=wsname, DataX=x, DataY=y, DataE=e, Nspec=nspec, UnitX=xunit, VerticalAxisUnit=vunit, VerticalAxisValues=vvalue) #add sample logs to new workspace - CopyLogs(InputWorkspace=temp_ws, OutputWorkspace=wsname) - addElwinLogs(wsname, label, eRange, range2_a) + CopyLogs(InputWorkspace=tempWS, OutputWorkspace=wsname) + addElwinLogs(wsname, label, eRange, Range2) # remove the temp workspace now we've copied the logs - DeleteWorkspace(temp_ws) + DeleteWorkspace(tempWS) if unit[0] == 'Temperature': - AddSampleLog(Workspace=e1_ws, LogName="temp_normalise", + AddSampleLog(Workspace=e1WS, LogName="temp_normalise", LogType="String", LogText=str(Normalise)) #create workspace normalized to the lowest temperature if Normalise: - elt_ws = ename+'_elt' - + eltWS = ename+'_elt' + #create elt workspace - mtd[elf_ws].clone(OutputWorkspace=elt_ws) - elwinNormalizeToLowestTemp(elt_ws) + mtd[elfWS].clone(OutputWorkspace=eltWS) + elwinNormalizeToLowestTemp(eltWS) #set labels and meta data - unitx = mtd[elt_ws].getAxis(0).setUnit("Label") + unitx = mtd[eltWS].getAxis(0).setUnit("Label") unitx.setLabel(unit[0], unit[1]) - addElwinLogs(elt_ws, label, eRange, range2_a) + addElwinLogs(eltWS, label, eRange, Range2) #append workspace name to output files list - wsnames.append(elt_ws) + wsnames.append(eltWS) #set labels on workspace axes - unity = mtd[e1_ws].getAxis(1).setUnit("Label") + unity = mtd[e1WS].getAxis(1).setUnit("Label") unity.setLabel(unit[0], unit[1]) - - unity = mtd[e2_ws].getAxis(1).setUnit("Label") + + unity = mtd[e2WS].getAxis(1).setUnit("Label") unity.setLabel(unit[0], unit[1]) - - unitx = mtd[elf_ws].getAxis(0).setUnit("Label") + + unitx = mtd[elfWS].getAxis(0).setUnit("Label") unitx.setLabel(unit[0], unit[1]) if Save: elwinSaveWorkspaces(wsnames, workdir, Verbose) if Plot: - elwinPlot(label, e1_ws, e2_ws, elf_ws, elt_ws) + elwinPlot(label,e1WS,e2WS,elfWS,eltWS) EndTime('Elwin') - return e1_ws, e2_ws + return e1WS,e2WS #normalize workspace to the lowest temperature def elwinNormalizeToLowestTemp(eltWS): - num_hist = mtd[eltWS].getNumberHistograms() - + nhist = mtd[eltWS].getNumberHistograms() + #normalize each spectrum in the workspace - for hist in range(0, num_hist): - y_vals = mtd[eltWS].readY(hist) - scale = 1.0 / y_vals[0] - y_vals_scaled = scale * y_vals - mtd[eltWS].setY(hist, y_vals_scaled) + for n in range(0,nhist): + y = mtd[eltWS].readY(n) + scale = 1.0/y[0] + yscaled = scale * y + mtd[eltWS].setY(n, yscaled) # Write each of the created workspaces to file def elwinSaveWorkspaces(flist, dir, Verbose): @@ -417,7 +418,7 @@ def addElwinLogs(ws, label, eRange, Range2): AddSampleLog(Workspace=ws, LogName="range2_end", LogType="Number", LogText=str(eRange[3])) #Plot each of the workspace output by elwin -def elwinPlot(label, eq1, eq2, elf, elt): +def elwinPlot(label,eq1,eq2,elf,elt): plotElwinWorkspace(eq1, yAxisTitle='Elastic Intensity', setScale=True) plotElwinWorkspace(eq2, yAxisTitle='log(Elastic Intensity)', setScale=True) plotElwinWorkspace(elf, xAxisTitle=label) @@ -427,67 +428,67 @@ def elwinPlot(label, eq1, eq2, elf, elt): #Plot a workspace generated by Elwin def plotElwinWorkspace(ws, xAxisTitle=None, yAxisTitle=None, setScale=False): - workspace = mtd[ws] - num_bins = ws.blocksize() - last_x = ws.readX(0)[num_bins-1] + ws = mtd[ws] + nBins = ws.blocksize() + lastX = ws.readX(0)[nBins-1] - num_hist = workspace.getNumberHistograms() + nhist = ws.getNumberHistograms() try: - graph = MTD_PLOT.plotSpectrum(ws, range(0, num_hist)) - except RuntimeError, _: + graph = mp.plotSpectrum(ws, range(0,nhist)) + except RuntimeError, e: #User clicked cancel on plot so don't do anything return None - + layer = graph.activeLayer() #set the x scale of the layer if setScale: - layer.setScale(MTD_PLOT.Layer.Bottom, 0.0, last_x) - + layer.setScale(mp.Layer.Bottom, 0.0, lastX) + #set the title on the x-axis if xAxisTitle: - layer.setAxisTitle(MTD_PLOT.Layer.Bottom, xAxisTitle) + layer.setAxisTitle(mp.Layer.Bottom, xAxisTitle) #set the title on the y-axis if yAxisTitle: - layer.setAxisTitle(MTD_PLOT.Layer.Left, yAxisTitle) + layer.setAxisTitle(mp.Layer.Left, yAxisTitle) ############################################################################## # Fury ############################################################################## def furyPlot(inWS, spec): - graph = MTD_PLOT.plotSpectrum(inWS, spec) + graph = mp.plotSpectrum(inWS, spec) layer = graph.activeLayer() - layer.setScale(MTD_PLOT.Layer.Left, 0, 1.0) + layer.setScale(mp.Layer.Left, 0, 1.0) def fury(samWorkspaces, res_workspace, rebinParam, RES=True, Save=False, Verbose=False, - Plot=False): - + Plot=False): + StartTime('Fury') workdir = config['defaultsave.directory'] - sam_temp = samWorkspaces[0] - _, npt = CheckHistZero(sam_temp) - x_in = mtd[sam_temp].readX(0) - d1 = x_in[1] - x_in[0] + samTemp = samWorkspaces[0] + nsam,npt = CheckHistZero(samTemp) + Xin = mtd[samTemp].readX(0) + d1 = Xin[1]-Xin[0] if d1 < 1e-8: error = 'Data energy bin is zero' logger.notice('ERROR *** ' + error) sys.exit(error) - d2 = x_in[npt-1] - x_in[npt-2] - dmin = min(d1, d2) + d2 = Xin[npt-1]-Xin[npt-2] + dmin = min(d1,d2) pars = rebinParam.split(',') - if float(pars[1]) <= dmin: + if (float(pars[1]) <= dmin): error = 'EWidth = ' + pars[1] + ' < smallest Eincr = ' + str(dmin) logger.notice('ERROR *** ' + error) sys.exit(error) - out_ws_list = [] + outWSlist = [] # Process RES Data Only Once - CheckAnalysers(sam_temp, res_workspace, Verbose) - nres, nptr = CheckHistZero(res_workspace) + CheckAnalysers(samTemp, res_workspace, Verbose) + nres,nptr = CheckHistZero(res_workspace) if nres > 1: - CheckHistSame(sam_temp, 'Sample', res_workspace, 'Resolution') + CheckHistSame(samTemp,'Sample', res_workspace, 'Resolution') tmp_res_workspace = '__tmp_' + res_workspace Rebin(InputWorkspace=res_workspace, OutputWorkspace=tmp_res_workspace, Params=rebinParam) @@ -495,14 +496,16 @@ def fury(samWorkspaces, res_workspace, rebinParam, RES=True, Save=False, Verbose ConvertToPointData(InputWorkspace=tmp_res_workspace, OutputWorkspace=tmp_res_workspace) ExtractFFTSpectrum(InputWorkspace=tmp_res_workspace, OutputWorkspace='res_fft', FFTPart=2) Divide(LHSWorkspace='res_fft', RHSWorkspace='res_int', OutputWorkspace='res') - for sam_ws in samWorkspaces: - Rebin(InputWorkspace=sam_ws, OutputWorkspace='sam_data', Params=rebinParam) + for samWs in samWorkspaces: + (direct, filename) = os.path.split(samWs) + (root, ext) = os.path.splitext(filename) + Rebin(InputWorkspace=samWs, OutputWorkspace='sam_data', Params=rebinParam) Integration(InputWorkspace='sam_data', OutputWorkspace='sam_int') ConvertToPointData(InputWorkspace='sam_data', OutputWorkspace='sam_data') ExtractFFTSpectrum(InputWorkspace='sam_data', OutputWorkspace='sam_fft', FFTPart=2) Divide(LHSWorkspace='sam_fft', RHSWorkspace='sam_int', OutputWorkspace='sam') # Create save file name - savefile = getWSprefix(sam_ws) + 'iqt' + savefile = getWSprefix(samWs) + 'iqt' outWSlist.append(savefile) Divide(LHSWorkspace='sam', RHSWorkspace='res', OutputWorkspace=savefile) #Cleanup Sample Files @@ -518,187 +521,190 @@ def fury(samWorkspaces, res_workspace, rebinParam, RES=True, Save=False, Verbose opath = os.path.join(workdir, savefile+'.nxs') # path name for nxs file SaveNexusProcessed(InputWorkspace=savefile, Filename=opath) if Verbose: - logger.notice('Output file : '+opath) + logger.notice('Output file : '+opath) # Clean Up RES files DeleteWorkspace(tmp_res_workspace) DeleteWorkspace('res_int') DeleteWorkspace('res_fft') DeleteWorkspace('res') if Plot: - specrange = range(0, mtd[out_ws_list[0]].getNumberHistograms()) - furyPlot(out_ws_list, specrange) + specrange = range(0,mtd[outWSlist[0]].getNumberHistograms()) + furyPlot(outWSlist, specrange) EndTime('Fury') - return out_ws_list + return outWSlist ############################################################################## # FuryFit ############################################################################## -def furyfitSeq(inputWS, func, ftype, startx, endx, intensities_constrained=False, Save=False, Plot='None', Verbose=False): - StartTime('FuryFit') - num_hist = mtd[inputWS].getNumberHistograms() - - # name stem for generated workspace - output_workspace = getWSprefix(inputWS) + 'fury_' + ftype + "0_to_" + str(num_hist-1) - - fit_type = ftype[:-2] - if Verbose: - logger.notice('Option: '+fit_type) - logger.notice(func) - - tmp_fit_workspace = "__furyfit_fit_ws" - CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=tmp_fit_workspace, XMin=startx, XMax=endx) - ConvertToHistogram(tmp_fit_workspace, OutputWorkspace=tmp_fit_workspace) - convertToElasticQ(tmp_fit_workspace) - # build input string for PlotPeakByLogValue - input_str = [tmp_fit_workspace + ',i%d' % i for i in range(0, num_hist)] - input_str = ';'.join(input_str) - - PlotPeakByLogValue(Input=input_str, OutputWorkspace=output_workspace, Function=func, - StartX=startx, EndX=endx, FitType='Sequential', CreateOutput=True) - - # remove unsused workspaces - DeleteWorkspace(output_workspace + '_NormalisedCovarianceMatrices') - DeleteWorkspace(output_workspace + '_Parameters') - - fit_group = output_workspace + '_Workspaces' - params_table = output_workspace + '_Parameters' - RenameWorkspace(output_workspace, OutputWorkspace=params_table) - - #create *_Result workspace - result_workspace = output_workspace + "_Result" - parameter_names = ['A0', 'Intensity', 'Tau', 'Beta'] - convertParametersToWorkspace(params_table, "axis-1", parameter_names, result_workspace) - - #set x units to be momentum transfer - axis = mtd[result_workspace].getAxis(0) - axis.setUnit("MomentumTransfer") - - #process generated workspaces - wsnames = mtd[fit_group].getNames() - for i, workspace in enumerate(wsnames): - output_ws = output_workspace + '_%d_Workspace' % i - RenameWorkspace(workspace, OutputWorkspace=output_ws) - - sample_logs = {'start_x': startx, 'end_x': endx, 'fit_type': ftype, - 'intensities_constrained': intensities_constrained, 'beta_constrained': False} - - CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fit_group) - CopyLogs(InputWorkspace=inputWS, OutputWorkspace=result_workspace) - - addSampleLogs(fit_group, sample_logs) - addSampleLogs(result_workspace, sample_logs) - - if Save: - save_workspaces = [result_workspace, fit_group] - furyFitSaveWorkspaces(save_workspaces, Verbose) - - if Plot != 'None': - furyfitPlotSeq(result_workspace, Plot) - - EndTime('FuryFit') - return result_workspace +def furyfitSeq(inputWS, func, ftype, startx, endx, intensities_constrained=False, Save=False, Plot='None', Verbose=False): + + StartTime('FuryFit') + nHist = mtd[inputWS].getNumberHistograms() + + #name stem for generated workspace + output_workspace = getWSprefix(inputWS) + 'fury_' + ftype + "0_to_" + str(nHist-1) + + fitType = ftype[:-2] + if Verbose: + logger.notice('Option: '+fitType) + logger.notice(func) + + tmp_fit_workspace = "__furyfit_fit_ws" + CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=tmp_fit_workspace, XMin=startx, XMax=endx) + ConvertToHistogram(tmp_fit_workspace, OutputWorkspace=tmp_fit_workspace) + convertToElasticQ(tmp_fit_workspace) + + #build input string for PlotPeakByLogValue + input_str = [tmp_fit_workspace + ',i%d' % i for i in range(0,nHist)] + input_str = ';'.join(input_str) + + PlotPeakByLogValue(Input=input_str, OutputWorkspace=output_workspace, Function=func, + StartX=startx, EndX=endx, FitType='Sequential', CreateOutput=True) + + #remove unsused workspaces + DeleteWorkspace(output_workspace + '_NormalisedCovarianceMatrices') + DeleteWorkspace(output_workspace + '_Parameters') + + fit_group = output_workspace + '_Workspaces' + params_table = output_workspace + '_Parameters' + RenameWorkspace(output_workspace, OutputWorkspace=params_table) + + #create *_Result workspace + result_workspace = output_workspace + "_Result" + parameter_names = ['A0', 'Intensity', 'Tau', 'Beta'] + convertParametersToWorkspace(params_table, "axis-1", parameter_names, result_workspace) + + #set x units to be momentum transfer + axis = mtd[result_workspace].getAxis(0) + axis.setUnit("MomentumTransfer") + + #process generated workspaces + wsnames = mtd[fit_group].getNames() + params = [startx, endx, fitType] + for i, ws in enumerate(wsnames): + output_ws = output_workspace + '_%d_Workspace' % i + RenameWorkspace(ws, OutputWorkspace=output_ws) + + sample_logs = {'start_x': startx, 'end_x': endx, 'fit_type': ftype, + 'intensities_constrained': intensities_constrained, 'beta_constrained': False} + + CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fit_group) + CopyLogs(InputWorkspace=inputWS, OutputWorkspace=result_workspace) + + addSampleLogs(fit_group, sample_logs) + addSampleLogs(result_workspace, sample_logs) + + if Save: + save_workspaces = [result_workspace, fit_group] + furyFitSaveWorkspaces(save_workspaces, Verbose) + + if Plot != 'None' : + furyfitPlotSeq(result_workspace, Plot) + + EndTime('FuryFit') + return result_workspace def furyfitMult(inputWS, function, ftype, startx, endx, intensities_constrained=False, Save=False, Plot='None', Verbose=False): - StartTime('FuryFit Multi') - - num_hist = mtd[inputWS].getNumberHistograms() - output_workspace = getWSprefix(inputWS) + 'fury_1Smult_s0_to_' + str(num_hist-1) - - option = ftype[:-2] - if Verbose: - logger.notice('Option: '+option) - logger.notice('Function: '+function) - - # prepare input workspace for fitting - tmp_fit_workspace = "__furyfit_fit_ws" - CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=tmp_fit_workspace, XMin=startx, XMax=endx) - ConvertToHistogram(tmp_fit_workspace, OutputWorkspace=tmp_fit_workspace) - convertToElasticQ(tmp_fit_workspace) - - # fit multi-domian functino to workspace - multi_domain_func, kwargs = createFuryMultiDomainFunction(function, tmp_fit_workspace) - Fit(Function=multi_domain_func, InputWorkspace=tmp_fit_workspace, WorkspaceIndex=0, - Output=output_workspace, CreateOutput=True, **kwargs) - - params_table = output_workspace + '_Parameters' - transposeFitParametersTable(params_table) - - # set first column of parameter table to be axis values - ax = mtd[tmp_fit_workspace].getAxis(1) - axis_values = ax.extractValues() - for i, value in enumerate(axis_values): - mtd[params_table].setCell('axis-1', i, value) - - # convert parameters to matrix workspace - result_workspace = output_workspace + "_Result" - parameter_names = ['A0', 'Intensity', 'Tau', 'Beta'] - convertParametersToWorkspace(params_table, "axis-1", parameter_names, result_workspace) - - # set x units to be momentum transfer - axis = mtd[result_workspace].getAxis(0) - axis.setUnit("MomentumTransfer") - - result_workspace = output_workspace + '_Result' - fit_group = output_workspace + '_Workspaces' - - sample_logs = {'start_x': startx, 'end_x': endx, 'fit_type': ftype, - 'intensities_constrained': intensities_constrained, 'beta_constrained': True} - - CopyLogs(InputWorkspace=inputWS, OutputWorkspace=result_workspace) - CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fit_group) - - addSampleLogs(result_workspace, sample_logs) - addSampleLogs(fit_group, sample_logs) - - DeleteWorkspace(tmp_fit_workspace) - - if Save: - save_workspaces = [result_workspace] - furyFitSaveWorkspaces(save_workspaces, Verbose) - - if Plot != 'None': - furyfitPlotSeq(result_workspace, Plot) - - EndTime('FuryFit Multi') - return result_workspace + StartTime('FuryFit Multi') + + nHist = mtd[inputWS].getNumberHistograms() + output_workspace = getWSprefix(inputWS) + 'fury_1Smult_s0_to_' + str(nHist-1) + + option = ftype[:-2] + if Verbose: + logger.notice('Option: '+option) + logger.notice('Function: '+function) + + #prepare input workspace for fitting + tmp_fit_workspace = "__furyfit_fit_ws" + CropWorkspace(InputWorkspace=inputWS, OutputWorkspace=tmp_fit_workspace, XMin=startx, XMax=endx) + ConvertToHistogram(tmp_fit_workspace, OutputWorkspace=tmp_fit_workspace) + convertToElasticQ(tmp_fit_workspace) + + #fit multi-domian functino to workspace + multi_domain_func, kwargs = createFuryMultiDomainFunction(function, tmp_fit_workspace) + Fit(Function=multi_domain_func, InputWorkspace=tmp_fit_workspace, WorkspaceIndex=0, + Output=output_workspace, CreateOutput=True, **kwargs) + + params_table = output_workspace + '_Parameters' + transposeFitParametersTable(params_table) + + #set first column of parameter table to be axis values + ax = mtd[tmp_fit_workspace].getAxis(1) + axis_values = ax.extractValues() + for i, value in enumerate(axis_values): + mtd[params_table].setCell('axis-1', i, value) + + #convert parameters to matrix workspace + result_workspace = output_workspace + "_Result" + parameter_names = ['A0', 'Intensity', 'Tau', 'Beta'] + convertParametersToWorkspace(params_table, "axis-1", parameter_names, result_workspace) + + #set x units to be momentum transfer + axis = mtd[result_workspace].getAxis(0) + axis.setUnit("MomentumTransfer") + + result_workspace = output_workspace + '_Result' + fit_group = output_workspace + '_Workspaces' + + sample_logs = {'start_x': startx, 'end_x': endx, 'fit_type': ftype, + 'intensities_constrained': intensities_constrained, 'beta_constrained': True} + + CopyLogs(InputWorkspace=inputWS, OutputWorkspace=result_workspace) + CopyLogs(InputWorkspace=inputWS, OutputWorkspace=fit_group) + + addSampleLogs(result_workspace, sample_logs) + addSampleLogs(fit_group, sample_logs) + + DeleteWorkspace(tmp_fit_workspace) + + if Save: + save_workspaces = [result_workspace] + furyFitSaveWorkspaces(save_workspaces, Verbose) + + if Plot != 'None': + furyfitPlotSeq(result_workspace, Plot) + + EndTime('FuryFit Multi') + return result_workspace def createFuryMultiDomainFunction(function, input_ws): - multi = 'composite=MultiDomainFunction,NumDeriv=1;' - comp = '(composite=CompositeFunction,$domains=i;' + function + ');' - - ties = [] - kwargs = {} - num_spectra = mtd[input_ws].getNumberHistograms() - for i in range(0, num_spectra): - multi += comp - kwargs['WorkspaceIndex_' + str(i)] = i - - if i > 0: - kwargs['InputWorkspace_' + str(i)] = input_ws - - # tie beta for every spectrum - tie = 'f%d.f1.Beta=f0.f1.Beta' % i - ties.append(tie) - - ties = ','.join(ties) - multi += 'ties=(' + ties + ')' - - return multi, kwargs + multi= 'composite=MultiDomainFunction,NumDeriv=1;' + comp = '(composite=CompositeFunction,$domains=i;' + function + ');' + + ties = [] + kwargs = {} + num_spectra = mtd[input_ws].getNumberHistograms() + for i in range(0, num_spectra): + multi += comp + kwargs['WorkspaceIndex_' + str(i)] = i + + if i > 0: + kwargs['InputWorkspace_' + str(i)] = input_ws + + #tie beta for every spectrum + tie = 'f%d.f1.Beta=f0.f1.Beta' % i + ties.append(tie) + + ties = ','.join(ties) + multi += 'ties=(' + ties + ')' + + return multi, kwargs def furyFitSaveWorkspaces(save_workspaces, Verbose): - workdir = getDefaultWorkingDirectory() - for workspace in save_workspaces: - # save workspace to default directory - fpath = os.path.join(workdir, workspace+'.nxs') - SaveNexusProcessed(InputWorkspace=workspace, Filename=fpath) + workdir = getDefaultWorkingDirectory() + for ws in save_workspaces: + #save workspace to default directory + fpath = os.path.join(workdir, ws+'.nxs') + SaveNexusProcessed(InputWorkspace=ws, Filename=fpath) - if Verbose: - logger.notice(workspace + ' output to file : '+fpath) + if Verbose: + logger.notice(ws + ' output to file : '+fpath) def furyfitPlotSeq(ws, plot): @@ -706,7 +712,7 @@ def furyfitPlotSeq(ws, plot): param_names = ['Intensity', 'Tau', 'Beta'] else: param_names = [plot] - + plotParameters(ws, *param_names) @@ -715,12 +721,12 @@ def furyfitPlotSeq(ws, plot): ############################################################################## def msdfitPlotSeq(inputWS, xlabel): - workspace = mtd[inputWS+'_A1'] - if len(workspace.readX(0)) > 1: - msd_plot = MTD_PLOT.plotSpectrum(inputWS+'_A1', 0, True) + ws = mtd[inputWS+'_A1'] + if len(ws.readX(0)) > 1: + msd_plot = mp.plotSpectrum(inputWS+'_A1',0,True) msd_layer = msd_plot.activeLayer() - msd_layer.setAxisTitle(MTD_PLOT.Layer.Bottom, xlabel) - msd_layer.setAxisTitle(MTD_PLOT.Layer.Left, '') + msd_layer.setAxisTitle(mp.Layer.Bottom,xlabel) + msd_layer.setAxisTitle(mp.Layer.Left,'') def msdfit(ws, startX, endX, spec_min=0, spec_max=None, Save=False, Verbose=False, Plot=True): StartTime('msdFit') @@ -740,38 +746,38 @@ def msdfit(ws, startX, endX, spec_min=0, spec_max=None, Save=False, Verbose=Fals xlabel = ws_run.getLogData('vert_axis').value mname = ws[:-4] - msd_ws = mname + '_msd' + msdWS = mname+'_msd' - # fit line to each of the spectra + #fit line to each of the spectra function = 'name=LinearBackground, A0=0, A1=0' - input_params = [ws+',i%d' % i for i in xrange(spec_min, spec_max+1)] + input_params = [ ws+',i%d' % i for i in xrange(spec_min, spec_max+1)] input_params = ';'.join(input_params) - PlotPeakByLogValue(Input=input_params, OutputWorkspace=msd_ws, Function=function, + PlotPeakByLogValue(Input=input_params, OutputWorkspace=msdWS, Function=function, StartX=startX, EndX=endX, FitType='Sequential', CreateOutput=True) - DeleteWorkspace(msd_ws + '_NormalisedCovarianceMatrices') - DeleteWorkspace(msd_ws + '_Parameters') - msd_parameters = msd_ws +'_Parameters' - RenameWorkspace(msd_ws, OutputWorkspace=msd_parameters) + DeleteWorkspace(msdWS + '_NormalisedCovarianceMatrices') + DeleteWorkspace(msdWS + '_Parameters') + msd_parameters = msdWS+'_Parameters' + RenameWorkspace(msdWS, OutputWorkspace=msd_parameters) params_table = mtd[msd_parameters] - - # msd value should be positive, but the fit output is negative + + #msd value should be positive, but the fit output is negative msd = params_table.column('A1') for i, value in enumerate(msd): params_table.setCell('A1', i, value * -1) - # create workspaces for each of the parameters + #create workspaces for each of the parameters group = [] - ws_name = msd_ws + '_A0' + ws_name = msdWS + '_A0' group.append(ws_name) ConvertTableToMatrixWorkspace(msd_parameters, OutputWorkspace=ws_name, ColumnX='axis-1', ColumnY='A0', ColumnE='A0_Err') xunit = mtd[ws_name].getAxis(0).setUnit('Label') xunit.setLabel('Temperature', 'K') - ws_name = msd_ws + '_A1' + ws_name = msdWS + '_A1' group.append(ws_name) ConvertTableToMatrixWorkspace(msd_parameters, OutputWorkspace=ws_name, ColumnX='axis-1', ColumnY='A1', ColumnE='A1_Err') @@ -781,39 +787,39 @@ def msdfit(ws, startX, endX, spec_min=0, spec_max=None, Save=False, Verbose=Fals xunit = mtd[ws_name].getAxis(0).setUnit('Label') xunit.setLabel('Temperature', 'K') - GroupWorkspaces(InputWorkspaces=','.join(group), OutputWorkspace=msd_ws) + GroupWorkspaces(InputWorkspaces=','.join(group),OutputWorkspace=msdWS) - # add sample logs to output workspace - fit_workspaces = msd_ws + '_Workspaces' - CopyLogs(InputWorkspace=ws, OutputWorkspace=msd_ws) - AddSampleLog(Workspace=msd_ws, LogName="start_x", LogType="Number", LogText=str(startX)) - AddSampleLog(Workspace=msd_ws, LogName="end_x", LogType="Number", LogText=str(endX)) - CopyLogs(InputWorkspace=msd_ws + '_A0', OutputWorkspace=fit_workspaces) + #add sample logs to output workspace + fit_workspaces = msdWS + '_Workspaces' + CopyLogs(InputWorkspace=ws, OutputWorkspace=msdWS) + AddSampleLog(Workspace=msdWS, LogName="start_x", LogType="Number", LogText=str(startX)) + AddSampleLog(Workspace=msdWS, LogName="end_x", LogType="Number", LogText=str(endX)) + CopyLogs(InputWorkspace=msdWS + '_A0', OutputWorkspace=fit_workspaces) if Plot: - msdfitPlotSeq(msd_ws, xlabel) + msdfitPlotSeq(msdWS, xlabel) if Save: - msd_path = os.path.join(workdir, msd_ws+'.nxs') # path name for nxs file - SaveNexusProcessed(InputWorkspace=msd_ws, Filename=msd_path, Title=msd_ws) + msd_path = os.path.join(workdir, msdWS+'.nxs') # path name for nxs file + SaveNexusProcessed(InputWorkspace=msdWS, Filename=msd_path, Title=msdWS) if Verbose: logger.notice('Output msd file : '+msd_path) EndTime('msdFit') return fit_workspaces -def plotInput(inputfiles, spectra=[]): - one_spectra = False +def plotInput(inputfiles,spectra=[]): + OneSpectra = False if len(spectra) != 2: spectra = [spectra[0], spectra[0]] - one_spectra = True + OneSpectra = True workspaces = [] - for in_file in inputfiles: - root = LoadNexus(Filename=in_file) - if not one_spectra: - GroupDetectors(root, root, DetectorList=range(spectra[0], spectra[1]+1)) + for file in inputfiles: + root = LoadNexus(Filename=file) + if not OneSpectra: + GroupDetectors(root, root, DetectorList=range(spectra[0],spectra[1]+1) ) workspaces.append(root) if len(workspaces) > 0: - graph = MTD_PLOT.plotSpectrum(workspaces, 0) + graph = mp.plotSpectrum(workspaces,0) graph.activeLayer().setTitle(", ".join(workspaces)) ############################################################################## @@ -823,21 +829,17 @@ def plotInput(inputfiles, spectra=[]): def CubicFit(inputWS, spec, Verbose=False): '''Uses the Mantid Fit Algorithm to fit a quadratic to the inputWS parameter. Returns a list containing the fitted parameter values.''' - function = 'name=Quadratic, A0=1, A1=0, A2=0' - Fit(Function=function, InputWorkspace=inputWS, WorkspaceIndex=spec, + fit = Fit(Function=function, InputWorkspace=inputWS, WorkspaceIndex=spec, CreateOutput=True, Output='Fit') table = mtd['Fit_Parameters'] - - a0 = table.cell(0, 1) - a1 = table.cell(1, 1) - a2 = table.cell(2, 1) - a_abs = [a0, a1, a2] - + A0 = table.cell(0,1) + A1 = table.cell(1,1) + A2 = table.cell(2,1) + Abs = [A0, A1, A2] if Verbose: - logger.notice('Group '+str(spec)+' of '+inputWS+' ; fit coefficients are : '+str(Abs)) - - return a_abs + logger.notice('Group '+str(spec)+' of '+inputWS+' ; fit coefficients are : '+str(Abs)) + return Abs def subractCanWorkspace(sample, can, output_name, rebin_can=False): '''Subtract the can workspace from the sample workspace. @@ -864,71 +866,71 @@ def applyCorrections(inputWS, canWS, corr, rebin_can=False, Verbose=False): '''Through the PolynomialCorrection algorithm, makes corrections to the input workspace based on the supplied correction values.''' # Corrections are applied in Lambda (Wavelength) - - efixed = getEfixed(inputWS) # Get efixed + + efixed = getEfixed(inputWS) # Get efixed ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='Wavelength', EMode='Indirect', EFixed=efixed) - name_stem = corr[:-4] + nameStem = corr[:-4] if canWS != '': - (_, can_run) = getInstrRun(canWS) - corrections = [name_stem+'_ass', name_stem+'_assc', name_stem+'_acsc', name_stem+'_acc'] - corrected_ws = sam_name +'Correct_'+ can_run + (instr, can_run) = getInstrRun(canWS) + corrections = [nameStem+'_ass', nameStem+'_assc', nameStem+'_acsc', nameStem+'_acc'] + CorrectedWS = sam_name +'Correct_'+ can_run ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='Wavelength', EMode='Indirect', EFixed=efixed) else: - corrections = [name_stem+'_ass'] - corrected_ws = sam_name +'Corrected' - num_hist = mtd[inputWS].getNumberHistograms() + corrections = [nameStem+'_ass'] + CorrectedWS = sam_name +'Corrected' + nHist = mtd[inputWS].getNumberHistograms() # Check that number of histograms in each corrections workspace matches # that of the input (sample) workspace - for workspace in corrections: - if mtd[workspace].getNumberHistograms() != num_hist: - raise ValueError('Mismatch: num of spectra in '+workspace+' and inputWS') + for ws in corrections: + if ( mtd[ws].getNumberHistograms() != nHist ): + raise ValueError('Mismatch: num of spectra in '+ws+' and inputWS') # Workspaces that hold intermediate results - corrected_sample_ws = '__csam' - corrected_can_ws = '__ccan' - for i in range(0, num_hist): # Loop through each spectra in the inputWS - ExtractSingleSpectrum(InputWorkspace=inputWS, OutputWorkspace=corrected_sample_ws, + CorrectedSampleWS = '__csam' + CorrectedCanWS = '__ccan' + for i in range(0, nHist): # Loop through each spectra in the inputWS + ExtractSingleSpectrum(InputWorkspace=inputWS, OutputWorkspace=CorrectedSampleWS, WorkspaceIndex=i) - if len(corrections) == 1: - ass = CubicFit(corrections[0], i, Verbose) - PolynomialCorrection(InputWorkspace=corrected_sample_ws, OutputWorkspace=corrected_sample_ws, - Coefficients=ass, Operation='Divide') - if i == 0: - CloneWorkspace(InputWorkspace=corrected_sample_ws, OutputWorkspace=corrected_ws) + if ( len(corrections) == 1 ): + Ass = CubicFit(corrections[0], i, Verbose) + PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS, + Coefficients=Ass, Operation='Divide') + if ( i == 0 ): + CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS) else: - ConjoinWorkspaces(InputWorkspace1=corrected_ws, InputWorkspace2=corrected_sample_ws) + ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS) else: - ExtractSingleSpectrum(InputWorkspace=canWS, OutputWorkspace=corrected_can_ws, + ExtractSingleSpectrum(InputWorkspace=canWS, OutputWorkspace=CorrectedCanWS, WorkspaceIndex=i) - acc = CubicFit(corrections[3], i, Verbose) - PolynomialCorrection(InputWorkspace=corrected_can_ws, OutputWorkspace=corrected_can_ws, - Coefficients=acc, Operation='Divide') - acsc = CubicFit(corrections[2], i, Verbose) - PolynomialCorrection(InputWorkspace=corrected_can_ws, OutputWorkspace=corrected_can_ws, - Coefficients=acsc, Operation='Multiply') - - subractCanWorkspace(corrected_sample_ws, corrected_can_ws, corrected_sample_ws, rebin_can=rebin_can) - - assc = CubicFit(corrections[1], i, Verbose) - PolynomialCorrection(InputWorkspace=corrected_sample_ws, OutputWorkspace=corrected_sample_ws, - Coefficients=assc, Operation='Divide') - if i == 0: - CloneWorkspace(InputWorkspace=corrected_sample_ws, OutputWorkspace=corrected_ws) + Acc = CubicFit(corrections[3], i, Verbose) + PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS, + Coefficients=Acc, Operation='Divide') + Acsc = CubicFit(corrections[2], i, Verbose) + PolynomialCorrection(InputWorkspace=CorrectedCanWS, OutputWorkspace=CorrectedCanWS, + Coefficients=Acsc, Operation='Multiply') + + subractCanWorkspace(CorrectedSampleWS, CorrectedCanWS, CorrectedSampleWS, rebin_can=rebin_can) + + Assc = CubicFit(corrections[1], i, Verbose) + PolynomialCorrection(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedSampleWS, + Coefficients=Assc, Operation='Divide') + if ( i == 0 ): + CloneWorkspace(InputWorkspace=CorrectedSampleWS, OutputWorkspace=CorrectedWS) else: - ConjoinWorkspaces(InputWorkspace1=corrected_ws, InputWorkspace2=corrected_sample_ws, + ConjoinWorkspaces(InputWorkspace1=CorrectedWS, InputWorkspace2=CorrectedSampleWS, CheckOverlapping=False) - + ConvertUnits(InputWorkspace=inputWS, OutputWorkspace=inputWS, Target='DeltaE', EMode='Indirect', EFixed=efixed) - ConvertUnits(InputWorkspace=corrected_ws, OutputWorkspace=corrected_ws, Target='DeltaE', + ConvertUnits(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS, Target='DeltaE', EMode='Indirect', EFixed=efixed) - ConvertSpectrumAxis(InputWorkspace=corrected_ws, OutputWorkspace=corrected_ws+'_rqw', + ConvertSpectrumAxis(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_rqw', Target='ElasticQ', EMode='Indirect', EFixed=efixed) - RenameWorkspace(InputWorkspace=corrected_ws, OutputWorkspace=corrected_ws+'_red') - + RenameWorkspace(InputWorkspace=CorrectedWS, OutputWorkspace=CorrectedWS+'_red') + if canWS != '': ConvertUnits(InputWorkspace=canWS, OutputWorkspace=canWS, Target='DeltaE', EMode='Indirect', EFixed=efixed) @@ -936,16 +938,16 @@ def applyCorrections(inputWS, canWS, corr, rebin_can=False, Verbose=False): DeleteWorkspace('Fit_NormalisedCovarianceMatrix') DeleteWorkspace('Fit_Parameters') DeleteWorkspace('Fit_Workspace') - return corrected_ws - + return CorrectedWS + def abscorFeeder(sample, container, geom, useCor, corrections, Verbose=False, RebinCan=False, ScaleOrNotToScale=False, factor=1, Save=False, PlotResult='None', PlotContrib=False): '''Load up the necessary files and then passes them into the main applyCorrections routine.''' StartTime('ApplyCorrections') workdir = config['defaultsave.directory'] - s_hist, _ = CheckHistZero(sample) - + s_hist,sxlen = CheckHistZero(sample) + diffraction_run = checkUnitIs(sample, 'dSpacing') sam_name = getWSprefix(sample) ext = '_red' @@ -962,7 +964,7 @@ def abscorFeeder(sample, container, geom, useCor, corrections, Verbose=False, Re if diffraction_run and not checkUnitIs(container, 'dSpacing'): raise ValueError("Sample and Can must both have the same units.") - (_, can_run) = getInstrRun(container) + (instr, can_run) = getInstrRun(container) scaled_container = "__apply_corr_scaled_container" if ScaleOrNotToScale: @@ -986,97 +988,99 @@ def abscorFeeder(sample, container, geom, useCor, corrections, Verbose=False, Re cor_result = applyCorrections(sample, container, corrections, RebinCan, Verbose) rws = mtd[cor_result + ext] + outNm = cor_result + '_Result_' if Save: - cred_path = os.path.join(workdir, cor_result + ext + '.nxs') + cred_path = os.path.join(workdir,cor_result + ext + '.nxs') SaveNexusProcessed(InputWorkspace=cor_result + ext, Filename=cred_path) if Verbose: logger.notice('Output file created : '+cred_path) + calc_plot = [cor_result + ext, sample] res_plot = cor_result+'_rqw' else: - if scaled_container == '': + if ( scaled_container == '' ): sys.exit('ERROR *** Invalid options - nothing to do!') else: - sub_result = sam_name + 'Subtract_' + can_run + sub_result = sam_name +'Subtract_'+ can_run if Verbose: logger.notice('Subtracting '+container+' from '+sample) subractCanWorkspace(sample, scaled_container, sub_result, rebin_can=RebinCan) if not diffraction_run: - ConvertSpectrumAxis(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_rqw', + ConvertSpectrumAxis(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_rqw', Target='ElasticQ', EMode='Indirect', EFixed=efixed) RenameWorkspace(InputWorkspace=sub_result, OutputWorkspace=sub_result+'_red') rws = mtd[sub_result+'_red'] - out_nm = sub_result + '_Result_' + outNm= sub_result + '_Result_' if Save: - sred_path = os.path.join(workdir, sub_result + ext + '.nxs') + sred_path = os.path.join(workdir,sub_result + ext + '.nxs') SaveNexusProcessed(InputWorkspace=sub_result + ext, Filename=sred_path) if Verbose: logger.notice('Output file created : ' + sred_path) - + if not diffraction_run: res_plot = sub_result + '_rqw' else: res_plot = sub_result + '_red' - - if PlotResult != 'None': + + if (PlotResult != 'None'): plotCorrResult(res_plot, PlotResult) - if scaled_container != '': + if ( scaled_container != '' ): sws = mtd[sample] cws = mtd[scaled_container] names = 'Sample,Can,Calc' - + x_unit = 'DeltaE' if diffraction_run: x_unit = 'dSpacing' - + for i in range(0, s_hist): # Loop through each spectra in the inputWS - data_x = np.array(sws.readX(i)) - data_y = np.array(sws.readY(i)) - data_e = np.array(sws.readE(i)) - data_x = np.append(data_x, np.array(cws.readX(i))) - data_y = np.append(data_y, np.array(cws.readY(i))) - data_e = np.append(data_e, np.array(cws.readE(i))) - data_x = np.append(data_x, np.array(rws.readX(i))) - data_y = np.append(data_y, np.array(rws.readY(i))) - data_e = np.append(data_e, np.array(rws.readE(i))) - fout = out_nm + str(i) - - CreateWorkspace(OutputWorkspace=fout, DataX=data_x, DataY=data_y, DataE=data_e, + dataX = np.array(sws.readX(i)) + dataY = np.array(sws.readY(i)) + dataE = np.array(sws.readE(i)) + dataX = np.append(dataX,np.array(cws.readX(i))) + dataY = np.append(dataY,np.array(cws.readY(i))) + dataE = np.append(dataE,np.array(cws.readE(i))) + dataX = np.append(dataX,np.array(rws.readX(i))) + dataY = np.append(dataY,np.array(rws.readY(i))) + dataE = np.append(dataE,np.array(rws.readE(i))) + fout = outNm + str(i) + + CreateWorkspace(OutputWorkspace=fout, DataX=dataX, DataY=dataY, DataE=dataE, Nspec=3, UnitX=x_unit, VerticalAxisUnit='Text', VerticalAxisValues=names) if i == 0: group = fout else: group += ',' + fout - GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=out_nm[:-1]) + GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=outNm[:-1]) if PlotContrib: - plotCorrContrib(out_nm+'0', [0, 1, 2]) + plotCorrContrib(outNm+'0',[0,1,2]) if Save: - res_path = os.path.join(workdir, out_nm[:-1]+'.nxs') - SaveNexusProcessed(InputWorkspace=out_nm[:-1], Filename=res_path) + res_path = os.path.join(workdir,outNm[:-1]+'.nxs') + SaveNexusProcessed(InputWorkspace=outNm[:-1],Filename=res_path) if Verbose: logger.notice('Output file created : '+res_path) DeleteWorkspace(cws) EndTime('ApplyCorrections') -def plotCorrResult(inWS, PlotResult): - n_hist = mtd[inWS].getNumberHistograms() - if PlotResult == 'Spectrum' or PlotResult == 'Both': - if n_hist >= 10: # only plot up to 10 hists - n_hist = 10 +def plotCorrResult(inWS,PlotResult): + nHist = mtd[inWS].getNumberHistograms() + if (PlotResult == 'Spectrum' or PlotResult == 'Both'): + if nHist >= 10: #only plot up to 10 hists + nHist = 10 plot_list = [] - for i in range(0, n_hist): + for i in range(0, nHist): plot_list.append(i) - res_plot = MTD_PLOT.plotSpectrum(inWS, plot_list) - if PlotResult == 'Contour' or PlotResult == 'Both': - if n_hist >= 5: # needs at least 5 hists for a contour - MTD_PLOT.importMatrixWorkspace(inWS).plotGraph2D() + res_plot=mp.plotSpectrum(inWS,plot_list) + if (PlotResult == 'Contour' or PlotResult == 'Both'): + if nHist >= 5: #needs at least 5 hists for a contour + mp.importMatrixWorkspace(inWS).plotGraph2D() -def plotCorrContrib(plot_list, n): - con_plot = MTD_PLOT.plotSpectrum(plot_list, n) +def plotCorrContrib(plot_list,n): + con_plot=mp.plotSpectrum(plot_list,n) From 70536d79205610d6e4daa4c1ef564e2211401a92 Mon Sep 17 00:00:00 2001 From: Dan Nixon Date: Fri, 15 Aug 2014 15:32:20 +0100 Subject: [PATCH 2/3] Revert "Fix failing system tests" This reverts commit 0df8c857511aeb71574d921f7dad130e795919df. --- Code/Mantid/scripts/Inelastic/IndirectAbsCor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py b/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py index a4671b746f59..7b929f876663 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py +++ b/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py @@ -195,7 +195,7 @@ def AbsRun(inputWS, geom, beam, num_can, size, density, sigs, siga, avar, verbos NSpec=ndet, UnitX='Wavelength', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) - CreateWorkspace(OutputWorkspace=assc_ws, DataX=data_x, DataY=data_a2, + CreateWorkspace(OutputWorkspace=ass_ws, DataX=data_x, DataY=data_a2, NSpec=ndet, UnitX='Wavelength', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) From 9c70fb818c83bc6beac66ff78ebe489f43d91375 Mon Sep 17 00:00:00 2001 From: Dan Nixon Date: Fri, 15 Aug 2014 15:32:22 +0100 Subject: [PATCH 3/3] Revert "More general refactoring" This reverts commit 4aa47109518d9ad9bacd4150357d83d134f1cfdf. --- .../Mantid/scripts/Inelastic/IndirectBayes.py | 2 +- .../scripts/Inelastic/IndirectCommon.py | 98 ++++----- .../scripts/Inelastic/IndirectJumpFit.py | 191 +++++++++--------- 3 files changed, 130 insertions(+), 161 deletions(-) diff --git a/Code/Mantid/scripts/Inelastic/IndirectBayes.py b/Code/Mantid/scripts/Inelastic/IndirectBayes.py index 68c98ba1dfe6..dced31c18e5a 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectBayes.py +++ b/Code/Mantid/scripts/Inelastic/IndirectBayes.py @@ -972,4 +972,4 @@ def ResNormPlot(inputWS,Plot): s_plot=mp.plotSpectrum(sWS,0,False) if (Plot == 'Fit' or Plot == 'All'): fWS = inputWS + '_ResNorm_Fit' - f_plot=mp.plotSpectrum(fWS,0,False) + f_plot=mp.plotSpectrum(fWS,0,False) \ No newline at end of file diff --git a/Code/Mantid/scripts/Inelastic/IndirectCommon.py b/Code/Mantid/scripts/Inelastic/IndirectCommon.py index ff2b3d64d7c4..39bf14535a5a 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectCommon.py +++ b/Code/Mantid/scripts/Inelastic/IndirectCommon.py @@ -4,26 +4,21 @@ from IndirectImport import import_mantidplot -import os.path -import math -import datetime -import re +import os.path, math, datetime, re import numpy as np import itertools def StartTime(prog): logger.notice('----------') - message = 'Program ' + prog + ' started @ ' + str(datetime.datetime.now()) + message = 'Program ' + prog +' started @ ' + str(datetime.datetime.now()) logger.notice(message) - def EndTime(prog): - message = 'Program ' + prog + ' ended @ ' + str(datetime.datetime.now()) + message = 'Program ' + prog +' ended @ ' + str(datetime.datetime.now()) logger.notice(message) logger.notice('----------') - def loadInst(instrument): workspace = '__empty_' + instrument if not mtd.doesExist(workspace): @@ -31,7 +26,6 @@ def loadInst(instrument): idf = idf_dir + instrument + '_Definition.xml' LoadEmptyInstrument(Filename=idf, OutputWorkspace=workspace) - def loadNexus(filename): ''' Loads a Nexus file into a workspace with the name based on the @@ -42,7 +36,6 @@ def loadNexus(filename): LoadNexus(Filename=filename, OutputWorkspace=name) return name - def getInstrRun(ws_name): ''' Get the instrument name and run number from a workspace. @@ -53,7 +46,7 @@ def getInstrRun(ws_name): workspace = mtd[ws_name] run_number = str(workspace.getRunNumber()) if run_number == '0': - # attempt to parse run number off of name + #attempt to parse run number off of name match = re.match(r'([a-zA-Z]+)([0-9]+)', ws_name) if match: run_number = match.group(2) @@ -66,7 +59,6 @@ def getInstrRun(ws_name): instrument = instrument.lower() return instrument, run_number - def getWSprefix(wsname): ''' Returns a string of the form '__' on which @@ -85,7 +77,7 @@ def getWSprefix(wsname): (instrument, run_number) = getInstrRun(wsname) if facility == 'ILL': - run_name = instrument + '_' + run_number + run_name = instrument + '_'+ run_number else: run_name = instrument + run_number @@ -103,12 +95,10 @@ def getWSprefix(wsname): return prefix - def getEfixed(workspace, det_index=0): inst = mtd[workspace].getInstrument() return inst.getNumberParameter("efixed-val")[0] - def checkUnitIs(ws, unit_id, axis_index=0): """ Check that the workspace has the correct units by comparing @@ -118,7 +108,6 @@ def checkUnitIs(ws, unit_id, axis_index=0): unit = axis.getUnit() return unit.unitID() == unit_id - # Get the default save directory and check it's valid def getDefaultWorkingDirectory(): workdir = config['defaultsave.directory'] @@ -128,7 +117,6 @@ def getDefaultWorkingDirectory(): return workdir - def createQaxis(inputWS): result = [] workspace = mtd[inputWS] @@ -141,7 +129,7 @@ def createQaxis(inputWS): efixed = getEfixed(inputWS, i) detector = workspace.getDetector(i) theta = detector.getTwoTheta(sample_pos, beam_pos) / 2 - lamda = math.sqrt(81.787 / efixed) + lamda = math.sqrt(81.787/efixed) q_value = 4 * math.pi * math.sin(theta) / lamda result.append(q_value) else: @@ -157,7 +145,6 @@ def createQaxis(inputWS): result.append(float(axis.label(i))) return result - def GetWSangles(in_ws): nhist = mtd[in_ws].getNumberHistograms() # get no. of histograms/groups source_pos = mtd[in_ws].getInstrument().getSource().getPos() @@ -166,44 +153,39 @@ def GetWSangles(in_ws): angles = [] # will be list of angles for index in range(0, nhist): detector = mtd[in_ws].getDetector(index) # get index - two_theta = detector.getTwoTheta(sample_pos, beam_pos) * 180.0 / math.pi # calc angle + two_theta = detector.getTwoTheta(sample_pos, beam_pos)*180.0/math.pi # calc angle angles.append(two_theta) # add angle return angles - def GetThetaQ(workspace): efixed = getEfixed(workspace) - wavelas = math.sqrt(81.787 / efixed) # elastic wavelength - k0_val = 4.0 * math.pi / wavelas + wavelas = math.sqrt(81.787/efixed) # elastic wavelength + k0_val = 4.0*math.pi/wavelas theta = np.array(GetWSangles(workspace)) q_val = k0_val * np.sin(0.5 * np.radians(theta)) return theta, q_val - def ExtractFloat(data_string): """ Extract float values from an ASCII string""" values = data_string.split() values = map(float, values) return values - def ExtractInt(data_string): """ Extract int values from an ASCII string""" values = data_string.split() values = map(int, values) return values - def PadArray(inarray, nfixed): # pad a list to specified size - npt = len(inarray) - padding = nfixed - npt - outarray = [] - outarray.extend(inarray) - outarray += [0] * padding - return outarray - + npt = len(inarray) + padding = nfixed-npt + outarray = [] + outarray.extend(inarray) + outarray += [0]*padding + return outarray def CheckAnalysers(in1_ws, in2_ws, verbose): ws1 = mtd[in1_ws] @@ -213,25 +195,23 @@ def CheckAnalysers(in1_ws, in2_ws, verbose): analyser2 = ws2.getInstrument().getStringParameter('analyser')[0] reflection2 = ws2.getInstrument().getStringParameter('reflection')[0] if analyser1 != analyser2: - raise ValueError('Workspace ' + in1_ws + ' and ' + in2_ws + ' have different analysers') + raise ValueError('Workspace '+in1_ws+' and '+in2_ws+' have different analysers') elif reflection1 != reflection2: - raise ValueError('Workspace ' + in1_ws + ' and ' + in2_ws + ' have different reflections') + raise ValueError('Workspace '+in1_ws+' and '+in2_ws+' have different reflections') else: if verbose: - logger.notice('Analyser is ' + analyser1 + reflection1) - + logger.notice('Analyser is '+analyser1+reflection1) def CheckHistZero(in_ws): nhist = mtd[in_ws].getNumberHistograms() # no. of hist/groups in WS if nhist == 0: - raise ValueError('Workspace ' + in_ws + ' has NO histograms') + raise ValueError('Workspace '+in_ws+' has NO histograms') x_in = mtd[in_ws].readX(0) ntc = len(x_in) - 1 # no. points from length of x array if ntc == 0: - raise ValueError('Workspace ' + in_ws + ' has NO points') + raise ValueError('Workspace '+in_ws+' has NO points') return nhist, ntc - def CheckHistSame(in1_ws, name1, in2_ws, name2): num_hist1 = mtd[in1_ws].getNumberHistograms() # no. of hist/groups in WS1 x_range1 = mtd[in1_ws].readX(0) @@ -248,34 +228,31 @@ def CheckHistSame(in1_ws, name1, in2_ws, name2): name1, in1_ws, x_len1, name2, in2_ws, x_len2) raise ValueError(error_str) - def CheckXrange(x_range, range_type): if not ((len(x_range) == 2) or (len(x_range) == 4)): raise ValueError(range_type + ' - Range must contain either 2 or 4 numbers') for lower, upper in zip(x_range[::2], x_range[1::2]): if math.fabs(lower) < 1e-5: - raise ValueError(range_type + ' - input minimum (' + str(lower) + ') is Zero') + raise ValueError(range_type+' - input minimum ('+str(lower)+') is Zero') if math.fabs(upper) < 1e-5: - raise ValueError(range_type + ' - input maximum (' + str(upper) + ') is Zero') + raise ValueError(range_type+' - input maximum ('+str(upper)+') is Zero') if upper < lower: - raise ValueError(range_type + ' - input max (' + str(upper) + ') < min (' + str(lower) + ')') - + raise ValueError(range_type+' - input max ('+str(upper)+') < min ('+str(lower)+')') def CheckElimits(e_range, x_in): x_len = len(x_in) - 1 if math.fabs(e_range[0]) < 1e-5: - raise ValueError('Elimits - input emin ( ' + str(e_range[0]) + ' ) is Zero') + raise ValueError('Elimits - input emin ( '+str(e_range[0])+' ) is Zero') if e_range[0] < x_in[0]: - raise ValueError('Elimits - input emin ( ' + str(e_range[0]) + ' ) < data emin ( ' + str(x_in[0]) + ' )') + raise ValueError('Elimits - input emin ( '+str(e_range[0])+' ) < data emin ( '+str(x_in[0])+' )') if math.fabs(e_range[1]) < 1e-5: - raise ValueError('Elimits - input emax ( ' + str(e_range[1]) + ' ) is Zero') + raise ValueError('Elimits - input emax ( '+str(e_range[1])+' ) is Zero') if e_range[1] > x_in[x_len]: - raise ValueError('Elimits - input emax ( ' + str(e_range[1]) + ' ) > data emax ( ' + str(x_in[x_len]) + ' )') + raise ValueError('Elimits - input emax ( '+str(e_range[1])+' ) > data emax ( '+str(x_in[x_len])+' )') if e_range[1] < e_range[0]: - raise ValueError('Elimits - input emax ( ' + str(e_range[1]) + ' ) < emin ( ' + str(e_range[0]) + ' )') - + raise ValueError('Elimits - input emax ( '+str(e_range[1])+' ) < emin ( '+str(e_range[0])+' )') def getInstrumentParameter(ws, param_name): """Get an named instrument parameter from a workspace. @@ -302,7 +279,6 @@ def getInstrumentParameter(ws, param_name): return param - def plotSpectra(ws, y_axis_title, indicies=[]): """ Plot a selection of spectra given a list of indicies @@ -325,7 +301,6 @@ def plotSpectra(ws, y_axis_title, indicies=[]): # User clicked cancel on plot so don't do anything return - def plotParameters(ws, *param_names): """ Plot a number of spectra given a list of parameter names @@ -344,7 +319,6 @@ def plotParameters(ws, *param_names): if len(indicies) > 0: plotSpectra(ws, name, indicies) - def convertToElasticQ(input_ws, output_ws=None): """ Helper function to convert the spectrum axis of a sample to ElasticQ. @@ -370,7 +344,6 @@ def convertToElasticQ(input_ws, output_ws=None): else: raise RuntimeError('Input workspace must have either spectra or numeric axis.') - def transposeFitParametersTable(params_table, output_table=None): """ Transpose the parameter table created from a multi domain Fit. @@ -386,7 +359,7 @@ def transposeFitParametersTable(params_table, output_table=None): table_ws = '__tmp_table_ws' table_ws = CreateEmptyTableWorkspace(OutputWorkspace=table_ws) - param_names = params_table.column(0)[:-1] # -1 to remove cost function + param_names = params_table.column(0)[:-1] #-1 to remove cost function param_values = params_table.column(1)[:-1] param_errors = params_table.column(2)[:-1] @@ -408,11 +381,11 @@ def transposeFitParametersTable(params_table, output_table=None): table_ws.addColumn('double', error_name) # output parameter values to table row - for i in xrange(0, params_table.rowCount() - 1, num_params): - row_values = param_values[i:(i + num_params)] - row_errors = param_errors[i:(i + num_params)] + for i in xrange(0, params_table.rowCount()-1, num_params): + row_values = param_values[i:i+num_params] + row_errors = param_errors[i:i+num_params] row = [value for pair in zip(row_values, row_errors) for value in pair] - row = [i / num_params] + row + row = [i/num_params] + row table_ws.addRow(row) if output_table is None: @@ -449,7 +422,7 @@ def convertParametersToWorkspace(params_table, x_column, param_names, output_nam workspace_names = [] for param_name in param_names: column_names = search_for_fit_params(param_name, params_table) - column_error_names = search_for_fit_params(param_name + '_Err', params_table) + column_error_names = search_for_fit_params(param_name+'_Err', params_table) param_workspaces = [] for name, error_name in zip(column_names, column_error_names): ConvertTableToMatrixWorkspace(params_table, x_column, name, error_name, OutputWorkspace=name) @@ -470,7 +443,7 @@ def convertParametersToWorkspace(params_table, x_column, param_names, output_nam temp_workspaces.append(temp_peak_ws) # join all peaks into a single workspace - # TODO: I'm not sure exactly what this is supposed to do + ##TODO: I'm not sure exactly what this is supposed to do temp_workspace = temp_workspaces[0] for temp_ws in temp_workspaces[1:]: ConjoinWorkspaces(temp_workspace, temp_peak_ws, False) @@ -485,7 +458,6 @@ def convertParametersToWorkspace(params_table, x_column, param_names, output_nam mtd[output_name].replaceAxis(1, axis) - def addSampleLogs(ws, sample_logs): """ Add a dictionary of logs to a workspace. diff --git a/Code/Mantid/scripts/Inelastic/IndirectJumpFit.py b/Code/Mantid/scripts/Inelastic/IndirectJumpFit.py index 059822ac479e..b4d43481583e 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectJumpFit.py +++ b/Code/Mantid/scripts/Inelastic/IndirectJumpFit.py @@ -2,104 +2,101 @@ from mantid import config, logger, mtd from IndirectCommon import * from IndirectImport import import_mantidplot -import sys -import os.path -MTD_PLOT = import_mantidplot() - +import sys, os.path +mp = import_mantidplot() # Jump programs -def JumpRun(samWS, jumpFunc, width, qmin, qmax, Verbose=False, Plot=False, Save=False): - StartTime('Jump fit : '+jumpFunc+' ; ') - - workdir = getDefaultWorkingDirectory() - - # select the width we wish to fit - spectrum_ws = "__" + samWS - ExtractSingleSpectrum(InputWorkspace=samWS, OutputWorkspace=spectrum_ws, WorkspaceIndex=width) - - # convert to HWHM - Scale(InputWorkspace=spectrum_ws, Factor=0.5, OutputWorkspace=spectrum_ws) - - # crop the workspace between the given ranges - if Verbose: - logger.notice('Cropping from Q= ' + str(qmin) + ' to ' + str(qmax)) - - # give the user some extra infromation if required - if Verbose: - in_gr = mtd[samWS].getRun() - try: - log = in_gr.getLogData('fit_program') - if log: - val = log.value - logger.notice('Fit program was : '+val) - except RuntimeError: - # if we couldn't find the fit program, just pass - pass - - logger.notice('Parameters in ' + samWS) - - x_vals = mtd[samWS].readX(0) - xmax = x_vals[-1] - - # select fit function to use - if jumpFunc == 'CE': - # Chudley-Elliott: HWHM=(1-sin*(Q*L)/(Q*L))/Tau - # for Q->0 W=Q^2*L^2/(6*Tau) - - tval = 1.0 / xmax - lval = 1.5 - func = 'name=ChudleyElliot, Tau=' + str(tval) + ', L=' + str(lval) - - elif jumpFunc == 'HallRoss': - # Hall-Ross: HWHM=(1-exp(-L*Q^2))/Tau - # for Q->0 W=A*Q^2*r - - tval = 1.0 / xmax - lval = 1.5 - func = 'name=HallRoss, Tau=' + str(tval) + ', L=' + str(lval) - - elif jumpFunc == 'Fick': - # Fick: HWHM=D*Q^2 - - y_vals = mtd[samWS].readY(0) - diff = (y_vals[2]-y_vals[0]) / ((x_vals[2]-x_vals[0]) * (x_vals[2]-x_vals[0])) - func = 'name=FickDiffusion, D='+str(diff) - - elif jumpFunc == 'Teixeira': - # Teixeira: HWHM=Q^2*L/((1+Q^2*L)*tau) - # for Q->0 W= - - tval = 1.0/xmax - lval = 1.5 - func = 'name=TeixeiraWater, Tau='+str(tval)+', L='+str(lval) - - else: - sys.exit("Error in Jump Fit: Invalid fit function supplied.") - return - - # run fit function - fit_workspace_base = samWS[:-10] + '_' + jumpFunc + 'fit' - Fit(Function=func, InputWorkspace=spectrum_ws, CreateOutput=True, Output=fit_workspace_base, StartX=qmin, EndX=qmax) - fit_workspace = fit_workspace_base + '_Workspace' - - CopyLogs(InputWorkspace=samWS, OutputWorkspace=fit_workspace) - AddSampleLog(Workspace=fit_workspace, LogName="jump_function", LogType="String", LogText=jumpFunc) - AddSampleLog(Workspace=fit_workspace, LogName="q_min", LogType="Number", LogText=str(qmin)) - AddSampleLog(Workspace=fit_workspace, LogName="q_max", LogType="Number", LogText=str(qmax)) - - # process output options - if Save: - fit_path = os.path.join(workdir, fit_workspace+'.nxs') - SaveNexusProcessed(InputWorkspace=fit_workspace, Filename=fit_path) - - if Verbose: - logger.notice('Fit file is ' + fit_path) - - if Plot: - JumpPlot(fit_workspace) - - EndTime('Jump fit : '+jumpFunc+' ; ') - +def JumpRun(samWS,jumpFunc,width,qmin,qmax,Verbose=False,Plot=False,Save=False): + StartTime('Jump fit : '+jumpFunc+' ; ') + + workdir = getDefaultWorkingDirectory() + + #select the width we wish to fit + spectrumWs = "__" + samWS + ExtractSingleSpectrum(InputWorkspace=samWS, OutputWorkspace=spectrumWs, WorkspaceIndex=width) + + #convert to HWHM + Scale(InputWorkspace=spectrumWs, Factor=0.5, OutputWorkspace=spectrumWs) + + #crop the workspace between the given ranges + if Verbose: + logger.notice('Cropping from Q= ' + str(qmin) +' to '+ str(qmax)) + + #give the user some extra infromation if required + if Verbose: + inGR = mtd[samWS].getRun() + try: + log = inGR.getLogData('fit_program') + if log: + val = log.value + logger.notice('Fit program was : '+val) + except RuntimeError: + #if we couldn't find the fit program, just pass + pass + + logger.notice('Parameters in ' + samWS) + + x = mtd[samWS].readX(0) + xmax = x[-1] + + #select fit function to use + if jumpFunc == 'CE': + # Chudley-Elliott: HWHM=(1-sin*(Q*L)/(Q*L))/Tau + # for Q->0 W=Q^2*L^2/(6*Tau) + + tval = 1.0/xmax + lval = 1.5 + func = 'name=ChudleyElliot, Tau='+str(tval)+', L='+str(lval) + + elif jumpFunc == 'HallRoss': + # Hall-Ross: HWHM=(1-exp(-L*Q^2))/Tau + # for Q->0 W=A*Q^2*r + + tval = 1.0/xmax + lval = 1.5 + func = 'name=HallRoss, Tau='+str(tval)+', L='+str(lval) + + elif jumpFunc == 'Fick': + # Fick: HWHM=D*Q^2 + + y = mtd[samWS].readY(0) + diff = (y[2]-y[0])/((x[2]-x[0])*(x[2]-x[0])) + func = 'name=FickDiffusion, D='+str(diff) + + elif jumpFunc == 'Teixeira': + # Teixeira: HWHM=Q^2*L/((1+Q^2*L)*tau) + # for Q->0 W= + + tval = 1.0/xmax + lval = 1.5 + func = 'name=TeixeiraWater, Tau='+str(tval)+', L='+str(lval) + + else: + sys.exit("Error in Jump Fit: Invalid fit function supplied.") + return + + #run fit function + fit_workspace_base = samWS[:-10] +'_'+ jumpFunc +'fit' + Fit(Function=func, InputWorkspace=spectrumWs, CreateOutput=True, Output=fit_workspace_base, StartX=qmin, EndX=qmax) + fit_workspace = fit_workspace_base + '_Workspace' + + CopyLogs(InputWorkspace=samWS, OutputWorkspace=fit_workspace) + AddSampleLog(Workspace=fit_workspace, LogName="jump_function", LogType="String", LogText=jumpFunc) + AddSampleLog(Workspace=fit_workspace, LogName="q_min", LogType="Number", LogText=str(qmin)) + AddSampleLog(Workspace=fit_workspace, LogName="q_max", LogType="Number", LogText=str(qmax)) + + #process output options + if Save: + fit_path = os.path.join(workdir,fit_workspace+'.nxs') + SaveNexusProcessed(InputWorkspace=fit_workspace, Filename=fit_path) + + if Verbose: + logger.notice('Fit file is ' + fit_path) + + if Plot: + JumpPlot(fit_workspace) + + EndTime('Jump fit : '+jumpFunc+' ; ') def JumpPlot(inputWS): - MTD_PLOT.plotSpectrum(inputWS, [0, 1, 2], True) + mp.plotSpectrum(inputWS,[0,1,2],True) \ No newline at end of file