From bbe8ad916cc8c3fb55c8d3d079905fec100602b1 Mon Sep 17 00:00:00 2001 From: Dan Nixon Date: Mon, 18 Aug 2014 09:42:02 +0100 Subject: [PATCH] Refactoring IndirectAbsCor.py Refs #9344 --- .../scripts/Inelastic/IndirectAbsCor.py | 339 +++++++++--------- 1 file changed, 172 insertions(+), 167 deletions(-) diff --git a/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py b/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py index 34388afef2b6..a11885dbde2e 100644 --- a/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py +++ b/Code/Mantid/scripts/Inelastic/IndirectAbsCor.py @@ -1,9 +1,11 @@ # IDA F2PY Absorption Corrections Wrapper ## Handle selection of .pyd files for absorption corrections -import platform, sys + +import sys from IndirectImport import * + if is_supported_f2py_platform(): - cylabs = import_f2py("cylabs") + CYLABS = import_f2py("cylabs") else: unsupported_message() @@ -11,109 +13,113 @@ from mantid.simpleapi import * from mantid import config, logger, mtd import math, os.path, numpy as np -mp = import_mantidplot() - -def WaveRange(inWS, efixed): -# create a list of 10 equi-spaced wavelengths spanning the input data - oWS = '__WaveRange' - ExtractSingleSpectrum(InputWorkspace=inWS, OutputWorkspace=oWS, WorkspaceIndex=0) - ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength', - EMode='Indirect', EFixed=efixed) - Xin = mtd[oWS].readX(0) - xmin = mtd[oWS].readX(0)[0] - xmax = mtd[oWS].readX(0)[len(Xin)-1] + +MTD_PLOT = import_mantidplot() + +PICONV = math.pi / 180.0 + + +def WaveRange(in_ws, e_fixed): + # create a list of 10 equi-spaced wavelengths spanning the input data + out_ws = '__WaveRange' + ExtractSingleSpectrum(InputWorkspace=in_ws, OutputWorkspace=out_ws, WorkspaceIndex=0) + ConvertUnits(InputWorkspace=out_ws, OutputWorkspace=out_ws, Target='Wavelength', + EMode='Indirect', EFixed=e_fixed) + x_in = mtd[out_ws].readX(0) + x_min = mtd[out_ws].readX(0)[0] + x_max = mtd[out_ws].readX(0)[len(x_in)-1] ebin = 0.5 - nw1 = int(xmin/ebin) - nw2 = int(xmax/ebin)+1 - w1 = nw1*ebin - w2 = nw2*ebin + nw1 = int(x_min/ebin) + nw2 = int(x_max/ebin) + 1 + w1 = nw1 * ebin + w2 = nw2 * ebin wave = [] nw = 10 ebin = (w2-w1)/(nw-1) - for l in range(0,nw): + for l in range(0, nw): wave.append(w1+l*ebin) - DeleteWorkspace(oWS) + DeleteWorkspace(out_ws) return wave -def CheckSize(size,geom,ncan,Verbose): +def CheckSize(size, geom, num_can, verbose): if geom == 'cyl': if (size[1] - size[0]) < 1e-4: - error = 'Sample outer radius not > inner radius' + error = 'Sample outer radius not > inner radius' logger.notice('ERROR *** '+error) sys.exit(error) else: - if Verbose: + if verbose: message = 'Sam : inner radius = '+str(size[0])+' ; outer radius = '+str(size[1]) logger.notice(message) if geom == 'flt': if size[0] < 1e-4: - error = 'Sample thickness is zero' + error = 'Sample thickness is zero' logger.notice('ERROR *** '+error) sys.exit(error) else: - if Verbose: + if verbose: logger.notice('Sam : thickness = '+str(size[0])) - if ncan == 2: + if num_can == 2: if geom == 'cyl': if (size[2] - size[1]) < 1e-4: - error = 'Can inner radius not > sample outer radius' + error = 'Can inner radius not > sample outer radius' logger.notice('ERROR *** '+error) sys.exit(error) else: - if Verbose: + if verbose: message = 'Can : inner radius = '+str(size[1])+' ; outer radius = '+str(size[2]) logger.notice(message) if geom == 'flt': if size[1] < 1e-4: - error = 'Can thickness is zero' + error = 'Can thickness is zero' logger.notice('ERROR *** '+error) sys.exit(error) else: - if Verbose: + if verbose: logger.notice('Can : thickness = '+str(size[1])) -def CheckDensity(density,ncan): +def CheckDensity(density, num_can): if density[0] < 1e-5: - error = 'Sample density is zero' + error = 'Sample density is zero' logger.notice('ERROR *** '+error) sys.exit(error) - if ncan == 2: + if num_can == 2: if density[1] < 1e-5: - error = 'Can density is zero' + error = 'Can density is zero' logger.notice('ERROR *** '+error) sys.exit(error) -def AbsRun(inputWS, geom, beam, ncan, size, density, sigs, siga, avar, Verbose, Save): +def AbsRun(inputWS, geom, beam, num_can, size, density, sigs, siga, avar, verbose, save): workdir = getDefaultWorkingDirectory() - if Verbose: + if verbose: logger.notice('Sample run : '+inputWS) # check that there is data - Xin = mtd[inputWS].readX(0) - if len(Xin) == 0: - error = 'Sample file has no data' + x_in = mtd[inputWS].readX(0) + if len(x_in) == 0: + error = 'Sample file has no data' logger.notice('ERROR *** '+error) sys.exit(error) - CheckSize(size,geom,ncan,Verbose) - CheckDensity(density,ncan) + CheckSize(size, geom, num_can, verbose) + CheckDensity(density, num_can) det = GetWSangles(inputWS) ndet = len(det) efixed = getEfixed(inputWS) - wavelas = math.sqrt(81.787/efixed) # elastic wavelength - waves = WaveRange(inputWS, efixed) # get wavelengths + wavelas = math.sqrt(81.787/efixed) # elastic wavelength + waves = WaveRange(inputWS, efixed) # get wavelengths nw = len(waves) run_name = getWSprefix(inputWS) - - if Verbose: + + if verbose: message = 'Sam : sigt = '+str(sigs[0])+' ; siga = '+str(siga[0])+' ; rho = '+str(density[0]) logger.notice(message) - if ncan == 2: + if num_can == 2: message = 'Can : sigt = '+str(sigs[1])+' ; siga = '+str(siga[1])+' ; rho = '+str(density[1]) logger.notice(message) @@ -124,112 +130,113 @@ def AbsRun(inputWS, geom, beam, ncan, size, density, sigs, siga, avar, Verbose, message = 'Detector angles : '+str(ndet)+' from '+str(det[0])+' to '+str(det[ndet-1]) logger.notice(message) - + name = run_name + geom wrk = workdir + run_name - wrk.ljust(120,' ') - - dataA1 = [] - dataA2 = [] - dataA3 = [] - dataA4 = [] + wrk.ljust(120, ' ') + + data_a1 = [] + data_a2 = [] + data_a3 = [] + data_a4 = [] - #initially set errors to zero - eZero = np.zeros(nw) + # initially set errors to zero + e_zero = np.zeros(nw) for n in range(ndet): - #geometry is flat + # geometry is flat if geom == 'flt': angles = [avar, det[n]] - (A1,A2,A3,A4) = FlatAbs(ncan, size, density, sigs, siga, angles, waves) + (a1, a2, a3, a4) = FlatAbs(num_can, size, density, sigs, siga, angles, waves) kill = 0 #geometry is a cylinder elif geom == 'cyl': astep = avar if (astep) < 1e-5: - error = 'Step size is zero' + error = 'Step size is zero' logger.notice('ERROR *** '+error) sys.exit(error) - + nstep = int((size[1] - size[0])/astep) if nstep < 20: - error = 'Number of steps ( '+str(nstep)+' ) should be >= 20' + error = 'Number of steps ( '+str(nstep)+' ) should be >= 20' logger.notice('ERROR *** '+error) sys.exit(error) angle = det[n] - kill, A1, A2, A3, A4 = cylabs.cylabs(astep, beam, ncan, size, + kill, a1, a2, a3, a4 = CYLABS.cylabs(astep, beam, num_can, size, density, sigs, siga, angle, wavelas, waves, n, wrk, 0) if kill == 0: - if Verbose: + if verbose: logger.notice('Detector '+str(n)+' at angle : '+str(det[n])+' * successful') - dataA1 = np.append(dataA1,A1) - dataA2 = np.append(dataA2,A2) - dataA3 = np.append(dataA3,A3) - dataA4 = np.append(dataA4,A4) + data_a1 = np.append(data_a1, a1) + data_a2 = np.append(data_a2, a2) + data_a3 = np.append(data_a3, a3) + data_a4 = np.append(data_a4, a4) else: error = 'Detector '+str(n)+' at angle : '+str(det[n])+' *** failed : Error code '+str(kill) logger.notice('ERROR *** '+error) sys.exit(error) - dataX = waves * ndet - qAxis = createQaxis(inputWS) + data_x = waves * ndet + q_axis = createQaxis(inputWS) # Create the output workspaces - assWS = name + '_ass' - asscWS = name + '_assc' - acscWS = name + '_acsc' - accWS = name + '_acc' + ass_ws = name + '_ass' + assc_ws = name + '_assc' + acsc_ws = name + '_acsc' + acc_ws = name + '_acc' fname = name +'_Abs' - CreateWorkspace(OutputWorkspace=assWS, DataX=dataX, DataY=dataA1, + CreateWorkspace(OutputWorkspace=ass_ws, DataX=data_x, DataY=data_a1, NSpec=ndet, UnitX='Wavelength', - VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) + VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) - CreateWorkspace(OutputWorkspace=asscWS, DataX=dataX, DataY=dataA2, + CreateWorkspace(OutputWorkspace=assc_ws, DataX=data_x, DataY=data_a2, NSpec=ndet, UnitX='Wavelength', - VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) + VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) - CreateWorkspace(OutputWorkspace=acscWS, DataX=dataX, DataY=dataA3, + CreateWorkspace(OutputWorkspace=acsc_ws, DataX=data_x, DataY=data_a3, NSpec=ndet, UnitX='Wavelength', - VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) + VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) - CreateWorkspace(OutputWorkspace=accWS, DataX=dataX, DataY=dataA4, + CreateWorkspace(OutputWorkspace=acc_ws, DataX=data_x, DataY=data_a4, NSpec=ndet, UnitX='Wavelength', - VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) + VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=q_axis) - group = assWS +','+ asscWS +','+ acscWS +','+ accWS - GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname) + group = ass_ws + ',' + assc_ws + ',' + acsc_ws + ',' + acc_ws + GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=fname) # save output to file if required - if Save: - opath = os.path.join(workdir,fname+'.nxs') + if save: + opath = os.path.join(workdir, fname+'.nxs') SaveNexusProcessed(InputWorkspace=fname, Filename=opath) - if Verbose: + if verbose: logger.notice('Output file created : '+opath) - if ncan > 1: - return [assWS, asscWS, acscWS, accWS] + if num_can > 1: + return [ass_ws, assc_ws, acsc_ws, acc_ws] else: - return [assWS] + return [ass_ws] -def plotAbs(workspaces, plotOpt): - if ( plotOpt == 'None' ): return +def plotAbs(workspaces, plot_opt): + if plot_opt == 'None': + return - if ( plotOpt == 'Wavelength' or plotOpt == 'Both' ): - graph = mp.plotSpectrum(workspaces, 0) + if plot_opt == 'Wavelength' or plot_opt == 'Both': + graph = MTD_PLOT.plotSpectrum(workspaces, 0) - if ( plotOpt == 'Angle' or plotOpt == 'Both' ): - graph = mp.plotTimeBin(workspaces, 0) - graph.activeLayer().setAxisTitle(mp.Layer.Bottom, 'Angle') + if plot_opt == 'Angle' or plot_opt == 'Both': + graph = MTD_PLOT.plotTimeBin(workspaces, 0) + graph.activeLayer().setAxisTitle(MTD_PLOT.Layer.Bottom, 'Angle') -def AbsRunFeeder(inputWS, canWS, geom, ncan, size, avar, density, beam_width=None, sampleFormula=None, canFormula=None, sigs=None, siga=None, - plotOpt='None', Verbose=False,Save=False): +def AbsRunFeeder(input_ws, can_ws, geom, ncan, size, avar, density, beam_width=None, sample_formula=None, can_formula=None, sigs=None, siga=None, + plotOpt='None', verbose=False, save=False): """ Handles the feeding of input and plotting of output for the F2PY absorption correction routine. @@ -253,14 +260,14 @@ def AbsRunFeeder(inputWS, canWS, geom, ncan, size, avar, density, beam_width=Non StartTime('CalculateCorrections') CheckDensity(density, ncan) - #attempt to find beam width if none given + # attempt to find beam width if none given if beam_width is None: - beam_width = getInstrumentParameter(inputWS, 'Workflow.beam-width') + beam_width = getInstrumentParameter(input_ws, 'Workflow.beam-width') beam_width = float(beam_width) - #attempt to find beam height from parameter file + # attempt to find beam height from parameter file try: - beam_height = getInstrumentParameter(inputWS, 'Workflow.beam-height') + beam_height = getInstrumentParameter(input_ws, 'Workflow.beam-height') beam_height = float(beam_height) except ValueError: # fall back on default value for beam height @@ -273,14 +280,14 @@ def AbsRunFeeder(inputWS, canWS, geom, ncan, size, avar, density, beam_width=Non # beam[7:8] hsdown,hsup bottom and top of scattered beam from sample b. beam = [beam_height, 0.5 * beam_width, -0.5 * beam_width, (beam_width / 2), -(beam_width / 2), 0.0, beam_height, 0.0, beam_height] - if sampleFormula is None and (sigs is None or siga is None): + if sample_formula is None and (sigs is None or siga is None): raise ValueError("Either a formula for the sample or values for the cross sections must be supplied.") #set sample material based on input or formula - if sampleFormula is not None: - SetSampleMaterial(InputWorkspace=inputWS, ChemicalFormula=sampleFormula, SampleNumberDensity=density[0]) + if sample_formula is not None: + SetSampleMaterial(InputWorkspace=input_ws, ChemicalFormula=sample_formula, SampleNumberDensity=density[0]) - sample = mtd[inputWS].sample() + sample = mtd[input_ws].sample() sam_mat = sample.getMaterial() # total scattering x-section @@ -288,11 +295,11 @@ def AbsRunFeeder(inputWS, canWS, geom, ncan, size, avar, density, beam_width=Non # absorption x-section siga[0] = sam_mat.absorbXSection() - if canFormula is not None and ncan == 2: + if can_formula is not None and ncan == 2: #set can material based on input or formula - SetSampleMaterial(InputWorkspace=canWS, ChemicalFormula=canFormula, SampleNumberDensity=density[1]) + SetSampleMaterial(InputWorkspace=can_ws, ChemicalFormula=can_formula, SampleNumberDensity=density[1]) - can_sample = mtd[canWS].sample() + can_sample = mtd[can_ws].sample() can_mat = can_sample.getMaterial() # total scattering x-section for can @@ -302,20 +309,19 @@ def AbsRunFeeder(inputWS, canWS, geom, ncan, size, avar, density, beam_width=Non siga[1] = can_mat.absorbXSection() siga[2] = can_mat.absorbXSection() - workspaces = AbsRun(inputWS, geom, beam, ncan, size, density, - sigs, siga, avar, Verbose, Save) + workspaces = AbsRun(input_ws, geom, beam, ncan, size, density, + sigs, siga, avar, verbose, save) EndTime('CalculateCorrections') plotAbs(workspaces, plotOpt) - def FlatAbs(ncan, thick, density, sigs, siga, angles, waves): - """ + """ FlatAbs - calculate flat plate absorption factors - + For more information See: - - MODES User Guide: http://www.isis.stfc.ac.uk/instruments/iris/data-analysis/modes-v3-user-guide-6962.pdf - - C J Carlile, Rutherford Laboratory report, RL-74-103 (1974) + - MODES User Guide: http://www.isis.stfc.ac.uk/instruments/iris/data-analysis/modes-v3-user-guide-6962.pdf + - C J Carlile, Rutherford Laboratory report, RL-74-103 (1974) @param sigs - list of scattering cross-sections @param siga - list of absorption cross-sections @@ -325,12 +331,11 @@ def FlatAbs(ncan, thick, density, sigs, siga, angles, waves): @param angles - list of angles @param waves - list of wavelengths """ - PICONV = math.pi/180. #can angle and detector angle tcan1, theta1 = angles - canAngle = tcan1*PICONV - theta = theta1*PICONV + can_angle = tcan1 * PICONV + theta = theta1 * PICONV # tsec is the angle the scattered beam makes with the normal to the sample surface. tsec = theta1-tcan1 @@ -343,68 +348,68 @@ def FlatAbs(ncan, thick, density, sigs, siga, angles, waves): acc = np.ones(nlam) # case where tsec is close to 90 degrees. CALCULATION IS UNRELIABLE - if (abs(abs(tsec)-90.0) < 1.0): + if abs(abs(tsec)-90.0) < 1.0: #default to 1 for everything return ass, assc, acsc, acc else: #sample & can scattering x-section - sampleScatt, canScatt = sigs[:2] + sample_scatt, can_scatt = sigs[:2] #sample & can absorption x-section - sampleAbs, canAbs = siga[:2] - #sample & can density - sampleDensity, canDensity = density[:2] + sample_abs, can_abs = siga[:2] + #sample & can density + sample_density, can_density = density[:2] #thickness of the sample and can - samThickness, canThickness1, canThickness2 = thick - - tsec = tsec*PICONV + sam_thickness, can_thickness1, can_thickness2 = thick + + tsec = tsec * PICONV - sec1 = 1./math.cos(canAngle) + sec1 = 1./math.cos(can_angle) sec2 = 1./math.cos(tsec) #list of wavelengths waves = np.array(waves) #sample cross section - sampleXSection = (sampleScatt + sampleAbs * waves /1.8) * sampleDensity + sample_x_section = (sample_scatt + sample_abs * waves /1.8) * sample_density #vector version of fact - vecFact = np.vectorize(Fact) - fs = vecFact(sampleXSection, samThickness, sec1, sec2) + vec_fact = np.vectorize(Fact) + fs = vec_fact(sample_x_section, sam_thickness, sec1, sec2) - sampleSec1, sampleSec2 = calcThicknessAtSec(sampleXSection, samThickness, [sec1, sec2]) + sample_sec1, sample_sec2 = calcThicknessAtSec(sample_x_section, sam_thickness, [sec1, sec2]) - if (sec2 < 0.): - ass = fs / samThickness + if sec2 < 0.0: + ass = fs / sam_thickness else: - ass= np.exp(-sampleSec2) * fs / samThickness + ass = np.exp(-sample_sec2) * fs / sam_thickness - useCan = (ncan > 1) - if useCan: + use_can = (ncan > 1) + if use_can: #calculate can cross section - canXSection = (canScatt + canAbs * waves /1.8) * canDensity - assc, acsc, acc = calcFlatAbsCan(ass, canXSection, canThickness1, canThickness2, sampleSec1, sampleSec2, [sec1, sec2]) + can_x_section = (can_scatt + can_abs * waves /1.8) * can_density + assc, acsc, acc = calcFlatAbsCan(ass, can_x_section, can_thickness1, can_thickness2, sample_sec1, sample_sec2, [sec1, sec2]) return ass, assc, acsc, acc -def Fact(xSection,thickness,sec1,sec2): - S = xSection*thickness*(sec1-sec2) - F = 1.0 - if (S == 0.): - F = thickness +def Fact(x_section, thickness, sec1, sec2): + s = x_section * thickness * (sec1 - sec2) + f = 1.0 + if s == 0.0: + f = thickness else: - S = (1-math.exp(-S))/S - F = thickness*S - return F + s = (1 - math.exp(-s)) / s + f = thickness * s + return f -def calcThicknessAtSec(xSection, thickness, sec): +def calcThicknessAtSec(x_section, thickness, sec): sec1, sec2 = sec - thickSec1 = xSection * thickness * sec1 - thickSec2 = xSection * thickness * sec2 + thick_sec1 = x_section * thickness * sec1 + thick_sec2 = x_section * thickness * sec2 - return thickSec1, thickSec2 + return thick_sec1, thick_sec2 -def calcFlatAbsCan(ass, canXSection, canThickness1, canThickness2, sampleSec1, sampleSec2, sec): +def calcFlatAbsCan(ass, can_x_section, can_thickness1, can_thickness2, sample_sec1, sample_sec2, sec): assc = np.ones(ass.size) acsc = np.ones(ass.size) acc = np.ones(ass.size) @@ -412,36 +417,36 @@ def calcFlatAbsCan(ass, canXSection, canThickness1, canThickness2, sampleSec1, s sec1, sec2 = sec #vector version of fact - vecFact = np.vectorize(Fact) - f1 = vecFact(canXSection,canThickness1,sec1,sec2) - f2 = vecFact(canXSection,canThickness2,sec1,sec2) + vec_fact = np.vectorize(Fact) + f1 = vec_fact(can_x_section, can_thickness1, sec1, sec2) + f2 = vec_fact(can_x_section, can_thickness2, sec1, sec2) - canThick1Sec1, canThick1Sec2 = calcThicknessAtSec(canXSection, canThickness1, sec) - canThick2Sec1, canThick2Sec2 = calcThicknessAtSec(canXSection, canThickness2, sec) + can_thick1_sec1, can_thick1_sec2 = calcThicknessAtSec(can_x_section, can_thickness1, sec) + can_thick2_sec1, can_thick2_sec2 = calcThicknessAtSec(can_x_section, can_thickness2, sec) - if (sec2 < 0.): - val = np.exp(-(canThick1Sec1-canThick1Sec2)) + if sec2 < 0.0: + val = np.exp(-(can_thick1_sec1 - can_thick1_sec2)) assc = ass * val acc1 = f1 acc2 = f2 * val acsc1 = acc1 - acsc2 = acc2 * np.exp(-(sampleSec1-sampleSec2)) + acsc2 = acc2 * np.exp(-(sample_sec1 - sample_sec2)) else: - val = np.exp(-(canThick1Sec1+canThick2Sec2)) + val = np.exp(-(can_thick1_sec1 + can_thick2_sec2)) assc = ass * val - acc1 = f1 * np.exp(-(canThick1Sec2+canThick2Sec2)) + acc1 = f1 * np.exp(-(can_thick1_sec2 + can_thick2_sec2)) acc2 = f2 * val - acsc1 = acc1 * np.exp(-sampleSec2) - acsc2 = acc2 * np.exp(-sampleSec1) + acsc1 = acc1 * np.exp(-sample_sec2) + acsc2 = acc2 * np.exp(-sample_sec1) - canThickness = canThickness1+canThickness2 + can_thickness = can_thickness1 + can_thickness2 - if(canThickness > 0.): - acc = (acc1+acc2)/canThickness - acsc = (acsc1+acsc2)/canThickness + if can_thickness > 0.0: + acc = (acc1 + acc2) / can_thickness + acsc = (acsc1 + acsc2) / can_thickness - return assc, acsc, acc \ No newline at end of file + return assc, acsc, acc