Ticket #7969: IndirectCommon.py

File IndirectCommon.py, 11.5 KB (added by Samuel Jackson, 7 years ago)

IndirectCommon with added function.

Line 
1from mantid.simpleapi import *
2from mantid import config, logger
3import sys, platform, os.path, math, datetime, re
4   
5def StartTime(prog):
6    logger.notice('----------')
7    message = 'Program ' + prog +' started @ ' + str(datetime.datetime.now())
8    logger.notice(message)
9
10def EndTime(prog):
11    message = 'Program ' + prog +' ended @ ' + str(datetime.datetime.now())
12    logger.notice(message)
13    logger.notice('----------')
14
15def loadInst(instrument):   
16    ws = '__empty_' + instrument
17    if not mtd.doesExist(ws):
18        idf_dir = config['instrumentDefinition.directory']
19        idf = idf_dir + instrument + '_Definition.xml'
20        LoadEmptyInstrument(Filename=idf, OutputWorkspace=ws)
21
22def loadNexus(filename):
23    '''Loads a Nexus file into a workspace with the name based on the
24    filename. Convenience function for not having to play around with paths
25    in every function.'''
26    name = os.path.splitext( os.path.split(filename)[1] )[0]
27    LoadNexus(Filename=filename, OutputWorkspace=name)
28    return name
29   
30def getInstrRun(file):
31    mo = re.match('([a-zA-Z]+)([0-9]+)',file)
32    instr_and_run = mo.group(0)          # instr name + run number
33    instr = mo.group(1)                  # instrument prefix
34    run = mo.group(2)                    # run number as string
35    return instr,run
36
37def getWSprefix(wsname,runfile=None):
38    '''Returns a string of the form '<ins><run>_<analyser><refl>_' on which
39    all of our other naming conventions are built.
40    The workspace is used to get the instrument parameters. If the runfile
41    string is given it is expected to be a string with instrument prefix
42    and run number. If it is empty then the workspace name is assumed to
43    contain this information
44    '''
45    if wsname == '':
46        return ''
47    if runfile is None:
48        runfile = wsname
49    ws = mtd[wsname]
50    facility = config['default.facility']
51    ws_run = ws.getRun()
52    if 'facility' in ws_run:
53        facility = ws_run.getLogData('facility').value
54    if facility == 'ILL':               
55        inst = ws.getInstrument().getName()
56        runNo = ws.getRun()['run_number'].value
57        run_name = inst + '_'+ runNo
58    else:
59        (instr, run) = getInstrRun(runfile)
60        run_name = instr + run
61    try:
62        analyser = ws.getInstrument().getStringParameter('analyser')[0]
63        reflection = ws.getInstrument().getStringParameter('reflection')[0]
64    except IndexError:
65        analyser = ''
66        reflection = ''
67    prefix = run_name + '_' + analyser + reflection + '_'
68    return prefix
69
70def getEfixed(workspace, detIndex=0):
71    inst = mtd[workspace].getInstrument()
72    return inst.getNumberParameter("efixed-val")[0]
73
74def getRunTitle(workspace):
75    ws = mtd[workspace]
76    title = ws.getRun()['run_title'].value.strip()
77    runNo = ws.getRun()['run_number'].value
78    inst = ws.getInstrument().getName()
79    ins = config.getFacility().instrument(ins).shortName().lower()
80    valid = "-_.() %s%s" % (string.ascii_letters, string.digits)
81    title = ''.join(ch for ch in title if ch in valid)
82    title = ins + runNo + '-' + title
83    return title
84
85def createQaxis(inputWS):
86    result = []
87    ws = mtd[inputWS]
88    nHist = ws.getNumberHistograms()
89    if ws.getAxis(1).isSpectra():
90        inst = ws.getInstrument()
91        samplePos = inst.getSample().getPos()
92        beamPos = samplePos - inst.getSource().getPos()
93        for i in range(0,nHist):
94            efixed = getEfixed(inputWS, i)
95            detector = ws.getDetector(i)
96            theta = detector.getTwoTheta(samplePos, beamPos) / 2
97            lamda = math.sqrt(81.787/efixed)
98            q = 4 * math.pi * math.sin(theta) / lamda
99            result.append(q)
100    else:
101        axis = ws.getAxis(1)
102        msg = 'Creating Axis based on Detector Q value: '
103        if not axis.isNumeric():
104            msg += 'Input workspace must have either spectra or numeric axis.'
105            logger.notice(msg)
106            sys.exit(msg)
107        if ( axis.getUnit().unitID() != 'MomentumTransfer' ):
108            msg += 'Input must have axis values of Q'
109            logger.notice(msg)
110            sys.exit(msg)
111        for i in range(0, nHist):
112            result.append(float(axis.label(i)))
113    return result
114
115def GetWSangles(inWS,verbose=False):
116    nhist = mtd[inWS].getNumberHistograms()                                             # get no. of histograms/groups
117    sourcePos = mtd[inWS].getInstrument().getSource().getPos()
118    samplePos = mtd[inWS].getInstrument().getSample().getPos() 
119    beamPos = samplePos - sourcePos
120    angles = []                                                                         # will be list of angles
121    for index in range(0, nhist):
122        detector = mtd[inWS].getDetector(index)                                 # get index
123        twoTheta = detector.getTwoTheta(samplePos, beamPos)*180.0/math.pi               # calc angle
124        angles.append(twoTheta)                                         # add angle
125    return angles
126
127def GetThetaQ(inWS):
128    nhist = mtd[inWS].getNumberHistograms()                                             # get no. of histograms/groups
129    efixed = getEfixed(inWS)
130    wavelas = math.sqrt(81.787/efixed)                                     # elastic wavelength
131    k0 = 4.0*math.pi/wavelas
132    d2r = math.pi/180.0
133    sourcePos = mtd[inWS].getInstrument().getSource().getPos()
134    samplePos = mtd[inWS].getInstrument().getSample().getPos() 
135    beamPos = samplePos - sourcePos
136    theta = []
137    Q = []
138    for index in range(0,nhist):
139        detector = mtd[inWS].getDetector(index)                                 # get index
140        twoTheta = detector.getTwoTheta(samplePos, beamPos)*180.0/math.pi               # calc angle
141        theta.append(twoTheta)                                          # add angle
142        Q.append(k0*math.sin(0.5*twoTheta*d2r))
143    return theta,Q
144
145def ExtractFloat(a):                              #extract values from line of ascii
146    extracted = []
147    elements = a.split()                                                        #split line on spaces
148    for n in elements:
149        extracted.append(float(n))
150    return extracted                                 #values as list
151
152def ExtractInt(a):                              #extract values from line of ascii
153    extracted = []
154    elements = a.split()                                                        #split line on spaces
155    for n in elements:
156        extracted.append(int(n))
157    return extracted                                 #values as list
158
159def PadArray(inarray,nfixed):                   #pad a list to specified size
160        npt=len(inarray)
161        padding = nfixed-npt
162        outarray=[]
163        outarray.extend(inarray)
164        outarray +=[0]*padding
165        return outarray
166
167def CheckAnalysers(in1WS,in2WS,Verbose):
168    ws1 = mtd[in1WS]
169    a1 = ws1.getInstrument().getStringParameter('analyser')[0]
170    r1 = ws1.getInstrument().getStringParameter('reflection')[0]
171    ws2 = mtd[in2WS]
172    a2 = ws2.getInstrument().getStringParameter('analyser')[0]
173    r2 = ws2.getInstrument().getStringParameter('reflection')[0]
174    if a1 != a2:
175        error = 'Workspace '+in1WS+' and '+in2WS+' have different analysers'
176        logger.notice('ERROR *** '+error)
177        sys.exit(error)
178    elif r1 != r2:
179        error = 'Workspace '+in1WS+' and '+in2WS+' have different reflections'
180        logger.notice('ERROR *** '+error)
181        sys.exit(error)
182    else:
183        if Verbose:
184            logger.notice('Analyser is '+a1+r1)
185
186def CheckHistZero(inWS):
187    nhist = mtd[inWS].getNumberHistograms()       # no. of hist/groups in WS
188    if nhist == 0:
189        error = 'Workspace '+inWS+' has NO histograms'                 
190        logger.notice('ERROR *** ' + error)
191        sys.exit(error)
192    Xin = mtd[inWS].readX(0)
193    ntc = len(Xin)-1                                            # no. points from length of x array
194    if ntc == 0:
195        error = 'Workspace '+inWS+' has NO points'                     
196        logger.notice('ERROR *** ' + error)
197        sys.exit(error)
198    return nhist,ntc
199
200def CheckHistSame(in1WS,name1,in2WS,name2):
201    nhist1 = mtd[in1WS].getNumberHistograms()       # no. of hist/groups in WS1
202    X1 = mtd[in1WS].readX(0)
203    xlen1 = len(X1)
204    nhist2 = mtd[in2WS].getNumberHistograms()       # no. of hist/groups in WS2
205    X2 = mtd[in2WS].readX(0)
206    xlen2 = len(X2)
207    if nhist1 != nhist2:                                # check that no. groups are the same
208        e1 = name1+' ('+in1WS+') histograms (' +str(nhist1) + ')'
209        e2 = name2+' ('+in2WS+') histograms (' +str(nhist2) + ')'
210        error = e1 + ' not = ' + e2
211        logger.notice('ERROR *** ' + error)
212        sys.exit(error)
213    elif xlen1 != xlen2:
214        e1 = name1+' ('+in1WS+') array length (' +str(xlen1) + ')'
215        e2 = name2+' ('+in2WS+') array length (' +str(xlen2) + ')'
216        error = e1 + ' not = ' + e2
217        logger.notice('ERROR *** ' + error)
218        sys.exit(error)
219
220def CheckXrange(xrange,type):
221    if  not ( ( len(xrange) == 2 ) or ( len(xrange) == 4 ) ):
222        error = type + ' - Range must contain either 2 or 4 numbers'
223        logger.notice(error)
224        sys.exit(error)
225    if math.fabs(xrange[0]) < 1e-5:
226        error = type + ' - input minimum ('+str(xrange[0])+') is Zero'                 
227        logger.notice('ERROR *** ' + error)
228        sys.exit(error)
229    if math.fabs(xrange[1]) < 1e-5:
230        error = type + ' - input maximum ('+str(xrange[1])+') is Zero'                 
231        logger.notice('ERROR *** ' + error)
232        sys.exit(error)
233    if xrange[1] < xrange[0]:
234        error = type + ' - input max ('+str(xrange[1])+') < min ('+xrange[0]+')'                       
235        logger.notice('ERROR *** ' + error)
236        sys.exit(error)
237    if len(xrange) >2:
238        if math.fabs(xrange[2]) < 1e-5:
239            error = type + '2 - input minimum ('+str(xrange[2])+') is Zero'                     
240            logger.notice('ERROR *** ' + error)
241            sys.exit(error)
242        if math.fabs(xrange[3]) < 1e-5:
243            error = type + '2 - input maximum ('+str(xrange[3])+') is Zero'                     
244            logger.notice('ERROR *** ' + error)
245            sys.exit(error)
246        if xrange[3] < xrange[2]:
247            error = type + '2 - input max ('+str(xrange[3])+') < min ('+xrange[2]+')'                   
248            logger.notice('ERROR *** ' + error)
249            sys.exit(error)
250
251def CheckElimits(erange,Xin):
252    nx = len(Xin)-1
253    if math.fabs(erange[0]) < 1e-5:
254        error = 'Elimits - input emin ( '+str(erange[0])+' ) is Zero'                   
255        logger.notice('ERROR *** ' + error)
256        sys.exit(error)
257    if erange[0] < Xin[0]:
258        error = 'Elimits - input emin ( '+str(erange[0])+' ) < data emin ( '+str(Xin[0])+' )'           
259        logger.notice('ERROR *** ' + error)
260        sys.exit(error)
261    if math.fabs(erange[1]) < 1e-5:
262        error = 'Elimits - input emax ( '+str(erange[1])+' ) is Zero'                   
263        logger.notice('ERROR *** ' + error)
264        sys.exit(error)
265    if erange[1] > Xin[nx]:
266        error = 'Elimits - input emax ( '+str(erange[1])+' ) > data emax ( '+str(Xin[nx])+' )'                 
267        logger.notice('ERROR *** ' + error)
268        sys.exit(error)
269    if erange[1] < erange[0]:
270        error = 'Elimits - input emax ( '+str(erange[1])+' ) < emin ( '+erange[0]+' )'                 
271        logger.notice('ERROR *** ' + error)
272        sys.exit(error)
273
274def CopyLogHeader(inWS,outWS):   #copies header data from inWS to outWS
275    inGR = mtd[inWS].getRun()
276    outGR = mtd[outWS].getRun()
277    run_header = inGR.getLogData('run_header').value
278    AddSampleLog(Workspace=outWS, LogName='run_header', LogType='String', LogText=str(run_header))
279    run_start = inGR.getLogData('run_start').value
280    AddSampleLog(Workspace=outWS, LogName='run_start', LogType='String', LogText=str(run_start))
281    run_title = inGR.getLogData('run_title').value
282    AddSampleLog(Workspace=outWS, LogName='run_title', LogType='String', LogText=str(run_title))
283    user_name = inGR.getLogData('user_name').value
284    AddSampleLog(Workspace=outWS, LogName='user_name', LogType='String', LogText=str(user_name))