Skip to content

Commit

Permalink
Add code to get Q and theta for spec axis in Q
Browse files Browse the repository at this point in the history
Refs #10454
  • Loading branch information
DanNixon committed Oct 31, 2014
1 parent d48fb15 commit 26a974e
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 12 deletions.
9 changes: 4 additions & 5 deletions Code/Mantid/scripts/Inelastic/IndirectBayes.py
Expand Up @@ -185,10 +185,9 @@ def QLRun(program,samWS,resWS,resnormWS,erange,nbins,Fit,wfile,Loop,Verbose,Plot
logger.notice('Sample is ' + samWS)
logger.notice('Resolution is ' + resWS)

if facility == 'ISIS':
CheckAnalysers(samWS,resWS,Verbose)
efix = getEfixed(samWS)
theta,Q = GetThetaQ(samWS)
CheckAnalysers(samWS,resWS,Verbose)
efix = getEfixed(samWS)
theta, Q = GetThetaQ(samWS)

nsam,ntc = CheckHistZero(samWS)

Expand Down Expand Up @@ -972,4 +971,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)
33 changes: 26 additions & 7 deletions Code/Mantid/scripts/Inelastic/IndirectCommon.py
Expand Up @@ -157,15 +157,34 @@ def GetWSangles(inWS):
return angles

def GetThetaQ(ws):
nhist = mtd[ws].getNumberHistograms() # get no. of histograms/groups
efixed = getEfixed(ws)
wavelas = math.sqrt(81.787/efixed) # elastic wavelength
k0 = 4.0*math.pi/wavelas
"""
Returns the theta and elastic Q for each spectrum in a given workspace.
@param ws Wotkspace to get theta and Q for
@returns A tuple containing a list of theta values and a list of Q values
"""

theta = np.array(GetWSangles(ws))
Q = k0 * np.sin(0.5 * np.radians(theta))
eFixed = getEfixed(ws)
wavelas = math.sqrt(81.787 / eFixed) # Elastic wavelength
k0 = 4.0 * math.pi / wavelas

axis = mtd[ws].getAxis(1)

# If axis is in spec number need to retrieve angles and calculate Q
if axis.isSpectra():
theta = np.array(GetWSangles(ws))
q = k0 * np.sin(0.5 * np.radians(theta))

# If axis is in Q need to calculate back to angles and just return axis values
else if axis.isNumeric() and axis.getUnit().unitID() == 'MomentumTransfer':
q = axis.extractValues()
theta = 2.0 * np.degrees(np.arcsin(q / k0))

# Out of options here
else:
raise RuntimeError('Cannot get theta and Q for workspace %s' % ws)

return theta, Q
return theta, q

def ExtractFloat(data_string):
""" Extract float values from an ASCII string"""
Expand Down

0 comments on commit 26a974e

Please sign in to comment.