mantid
mantid copied to clipboard
mantidproject
Qsample on peak doesn't match actual Q
Expected behavior
That Qsample = 2 * pi * UB * HKL
Actual behavior
it looks like d = 2*pi/|Qsample| is only good to ~2 decimal points
Steps to reproduce the behavior
(1) generate peaks with a known easy to interpret UB and inspect Qsample in peak table vs calculated
# generate peak table
ws = LoadEmptyInstrument(InstrumentName='SXD', OutputWorkspace='ws')
axis = ws.getAxis(0)
axis.setUnit("TOF") # need this to add peak to table
UB = 0.5*np.eye(3)
SetUB(ws, UB=UB)
peaks = PredictPeaks(InputWorkspace=ws, OutputWorkspace='peaks')
mod = lambda v: np.sqrt(np.sum(np.array(v)**2))
norm = lambda v: np.array(v) / mod(v)
pk = peaks.getPeak(0)
qsam = pk.getQSampleFrame()
qcalc = 2*np.pi*peaks.sample().getOrientedLattice().getUB() @ pk.getHKL()
print('Qsample', qsam)
print("UB* HKL", qcalc)
print('Qsample unit vec', norm(qsam))
print("UB* HKL unit vec", norm(qcalc))
# Qsample [-0.0119224, 3.14155, 3.13119]
# UB* HKL [0, 3.14159265, 3.14159265] <- this is correct
# Qsample unit vec [-0.00268795, 0.70827101, 0.70593551]
# UB* HKL unit vec [0, 0.70710678, 0.70710678]
There should be 0 momentum transfer along x (i.e. first element = 0) and the component along y and z should be equal. The calculated q and the Qsample in the peak table should match but they don't. The
(2) Compare the mod Q
# UB correpsonds to unit cell
u = UnitCell(2,2,2,90,90,90) # corresponds to UB set
print(u.d(*pk.getHKL())) # 1.41421 <- this is the correct d spacing (i.e. sqrt(2))
print(2*np.pi/mod(qcalc))# 1.41421
print(pk.getDSpacing()) # 1.41655 <- this is what is in the peak table
print(2*np.pi/mod(qsam)) # 1.41655
(3) Look at d-spacing of all peaks - not consistent
Platforms affected
Is this because it works out Q by putting peak on nearest detector? I would still expect d-spacing to be the same as nominal.
Seems to affect other instruments, as well
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
# generate peak table
ws = LoadEmptyInstrument(InstrumentName='TOPAZ', OutputWorkspace='ws')
axis = ws.getAxis(0)
axis.setUnit("TOF") # need this to add peak to table
UB = 0.25*np.eye(3)
SetUB(ws, UB=UB)
peaks = PredictPeaks(InputWorkspace=ws, OutputWorkspace='peaks')
mod = lambda v: np.sqrt(np.sum(np.array(v)**2))
norm = lambda v: np.array(v) / mod(v)
pk = peaks.getPeak(0)
qsam = pk.getQSampleFrame()
qcalc = 2*np.pi*peaks.sample().getOrientedLattice().getUB() @ pk.getHKL()
print('Qsample', qsam)
print("UB* HKL", qcalc)
print('Qsample unit vec', norm(qsam))
print("UB* HKL unit vec", norm(qcalc))
# Qsample [-3.14086,1.57043,4.71514]
# UB* HKL [-3.14159265 1.57079633 4.71238898]
# Qsample unit vec [-0.5342413 0.26712006 0.80201815]
# UB* HKL unit vec [-0.53452248 0.26726124 0.80178373]
As observed with CORELLI data, the finite pixel size (and tube-gap) affect the predicted Q-lab vectors as the peak has to be tied to a pixel. The pixel coordinates is ultimately used to calculate the Q-vectors even when using PredictPeaks
Actually, trying your script with LeanElasticPeaks that doesn't have the pixel information, the result also doesn't agree.
# import mantid algorithms, numpy and matplotlib
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
# generate peak table
ws = LoadEmptyInstrument(InstrumentName='SXD', OutputWorkspace='ws')
axis = ws.getAxis(0)
axis.setUnit("TOF") # need this to add peak to table
UB = 0.5*np.eye(3)
SetUB(ws, UB=UB)
peaks = PredictPeaks(InputWorkspace=ws, OutputWorkspace='peaks', OutputType='LeanElasticPeak')
mod = lambda v: np.sqrt(np.sum(np.array(v)**2))
norm = lambda v: np.array(v) / mod(v)
pk = peaks.getPeak(2)
qsam = pk.getQSampleFrame()
qcalc = 2*np.pi*peaks.sample().getOrientedLattice().getUB() @ pk.getHKL()
print('Qsample', qsam)
print("UB* HKL", qcalc)
print('Qsample unit vec', norm(qsam))
print("UB* HKL unit vec", norm(qcalc))
# Qsample [0.0169505,3.1415,3.15854]
# UB* HKL [ 0. 3.14159265 3.14159265]
# Qsample unit vec [ 0.00380496 0.70518633 0.70901182]
# UB* HKL unit vec [ 0. 0.70710678 0.70710678]
Moved the milestone to 6.2 S3 so I don't forget about this but I don't think I'll find time to look into it. @zjmorgan would you have time to look into this?
@RichardWaiteSTFC I will look into this
That's great thanks! I don't know if it helps but from the look of this code we seem to calculate the d-spacing using the two-theta of the nominal detector (which in general will be slightly different from the two-theta of the peak, say if we've used predicted peaks rather than found peaks). For predicted peaks the d-spacing should not be calculated in this way. Maybe the best thing to do is when a peak is created (using DetectorID & TOF or QSample etc.) we should calculate d-spacing as a member variable in the creator and the getter should just retrieve it? https://github.com/mantidproject/mantid/blob/bb269d6b76dd029d3419f5df77cb73da6ccf2496/Framework/DataObjects/src/Peak.cpp#L421
Also QLab (from which QSample is retrieved) is also bound to the detector direction at least for Peaks https://github.com/mantidproject/mantid/blob/bb269d6b76dd029d3419f5df77cb73da6ccf2496/Framework/DataObjects/src/Peak.cpp#L446 This shouldn't be an issue for LeanPeaks though right?
@RichardWaiteSTFC That is correct. So the fact that it also affects the LeanPeaks suggests it might have some other cause within PredictPeaks.
@RichardWaiteSTFC That is correct. So the fact that it also affects the LeanPeaks suggests it might have some other cause within PredictPeaks.
Agreed, but we still should not be calculating d-spacing and Qlab for Peaks using the detector position as we know that's going to be wrong for predicted peaks...but as you say it looks like another issue
Also just found that when you create a peak with QSample (or equivalently HKL) it seems to move it to have QSample of nearest detector...
pk = peaks.createPeakQSample([0., 3.14159265, 3.14159265]) # same as first peak in peaks
print(pk.getQSampleFrame()) # [0.0169505, 3.1415, 3.15854]
pk = peaks.createPeakHKL([0,1,1]) # same as first peak in peaks
print(pk.getQSampleFrame()) # # [0.0169505, 3.1415, 3.15854]
Perhaps this is why the LeanPeaks are also incorrect - see this code (where q is QSample) https://github.com/mantidproject/mantid/blob/398fa38e84407d7c8f4f6787685bdd39e4cfcd7a/Framework/Crystal/src/PredictPeaks.cpp#L638 Don't know where the detector info would creep in, maybe it calls a function used for normal peaks at some point?
For normal peaks it uses the detector ID (which presumably is only the nearest detector to kf) and the wavelength https://github.com/mantidproject/mantid/blob/398fa38e84407d7c8f4f6787685bdd39e4cfcd7a/Framework/Crystal/src/PredictPeaks.cpp#L581 But there is some code to handle peaks that don't hit detectors (i.e. in extended detector space) - so perhaps that should be used instead even if it does hit a detector (that might screw other stuff up though)?
This really does seem to be the issue with nearest neighbor detector being used to calculate Qlab and Qsample. The LeanElasticWorkspace has this same behavior if setting CalculateWavelength=True. This will try to grab the pixel information to calculate the wavelength. If we don't care about the wavelength or energy (typically only true for monochromatic beam), we can set CalculateWavelength=False which will give the exact Q expected.
from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
# generate peak table
ws = LoadEmptyInstrument(InstrumentName='SXD', OutputWorkspace='ws')
axis = ws.getAxis(0)
axis.setUnit("TOF") # need this to add peak to table
UB = 0.5*np.eye(3)
SetUB(ws, UB=UB)
peaks = PredictPeaks(InputWorkspace=ws,
OutputWorkspace='peaks',
OutputType='LeanElasticPeak',
CalculateWavelength=False)
mod = lambda v: np.sqrt(np.sum(np.array(v)**2))
norm = lambda v: np.array(v) / mod(v)
pk = peaks.getPeak(2)
qsam = pk.getQSampleFrame()
qcalc = 2*np.pi*peaks.sample().getOrientedLattice().getUB() @ pk.getHKL()
print('Qsample', qsam)
print("UB* HKL", qcalc)
print('Qsample unit vec', norm(qsam))
print("UB* HKL unit vec", norm(qcalc))