Reduction#
The main data tasks with single crystal neutron diffraction involve loading data, determining and refining the UB-matrix, integrating peaks, and normalizing data into reciprocal space volumes.
Loading data#
The data files (nexus/*.nxs.h5) produced by each instrument using EPICS are similar.
Each of these files (typically corresponding to a fixed goniometer angle) can always be loaded with LoadEventNexus.
By loading the data, an EventWorkspace is created with time-of-flight (TOF) vs counts.
It is possible to ConvertUnits from time-of-flight (TOF) to some other physical unit.
Some examples include:
d-Spacing
Momentum
Momentum transfer
Wavelength
Of course, for monochromatic instruments like WAND², TOF is not especially useful as the incident wavelength is fixed.
Additionally, to collect a volume of reciprocal space, it is necessary to rotate the sample through an angular range.
This contrasts with Laue instruments where the wavelength band gives rise to a volume of reciprocal space for each angle.
Therefore, many files must be combined of detector counts vs rotation angle as an MDHistoWorkspace.
For WAND², the custom loader LoadWANDSCD can be used.
For DEMAND which uses SPICE, the experiment already produce a scan file with detector counts varying with rotation angle.
As such, files (.nxs) are auto-reduced and can be loaded with HB3AAdjustSampleNorm.
WAND²#
A salt crystal example (NaCl) described by space group Fm-3m.
### experiment data info ###
IPTS = 7776
runs = range(26640, 27944)
### -------------------- ###
### sample info ###
a = 5.64
b = 5.64
c = 5.64
alpha = 90
beta = 90
gamma = 90
centering = 'F'
### ----------- ###
### incident wavelength ###
wavelength = 1.488
### ------------------- ###
filename = '/HFIR/HB2C/IPTS-{}/nexus/HB2C_{}.nxs.h5'
filenames = ','.join([filename.format(IPTS, run) for run in runs])
LoadWANDSCD(Filename=filenames,
Grouping='4x4',
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='s1,0,1,0,1',
Average=False)
DEMAND#
The example is triphylite (LiFePO₄) with space group Pmnb.
### experiment data info ###
IPTS = 9884
exp = 817
scan = 2
### -------------------- ###
### sample info ###
a = 6.02
b = 10.35
c = 4.70
alpha = 90
beta = 90
gamma = 90
centering = 'P'
### ----------- ###
### incident wavelength ###
wavelength = 1.003
### ------------------- ###
wavelength_band = [0.98*wavelength, 1.02*wavelength]
filename = '/HFIR/HB3A/IPTS-{}/shared/autoreduce/HB3A_exp{:04}_scan{:04}.nxs'.format(IPTS, exp, scan)
HB3AAdjustSampleNorm(Filename=filename,
OutputType='Detector',
NormaliseBy='None',
Grouping='4x4',
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='omega,0,1,0,-1',
Axis1='chi,0,0,1,-1',
Axis2='phi,0,1,0,-1',
Average=False)
MANDI#
A syntentic garnet crystal (Yb₃Al₅O₁₂) that has space group Ia-3d.
### experiment data info ###
IPTS = 8776
run = 10934
### -------------------- ###
### sample info ###
a = 11.93
b = 11.93
c = 11.93
alpha = 90
beta = 90
gamma = 90
centering = 'I'
### ----------- ###
### incident wavelength ###
wavelength_band = [1, 3]
### ------------------- ###
filename = '/SNS/MANDI/IPTS-{}/nexus/MANDI_{}.nxs.h5'.format(IPTS, run)
LoadEventNexus(Filename=filename,
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='BL11B:Mot:omega,0,1,0,1',
Axis1='BL11B:Mot:chi,0,0,1,1',
Axis2='BL11B:Mot:phi,0,1,0,1',
Average=True)
PreprocessDetectorsToMD(InputWorkspace='data',
OutputWorkspace='detectors')
TOPAZ#
The same garnet calibration standard.
### experiment data info ###
IPTS = 31189
run = 46917
### -------------------- ###
### sample info ###
a = 11.93
b = 11.93
c = 11.93
alpha = 90
beta = 90
gamma = 90
centering = 'I'
### ----------- ###
### incident wavelength ###
wavelength_band = [0.3, 3.5]
### ------------------- ###
filename = '/SNS/TOPAZ/IPTS-{}/nexus/TOPAZ_{}.nxs.h5'.format(IPTS, run)
LoadEventNexus(Filename=filename,
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='BL12:Mot:omega,0,1,0,1',
Axis1='BL12:Mot:chi,0,0,1,1',
Axis2='BL12:Mot:phi,0,1,0,1',
Average=True)
PreprocessDetectorsToMD(InputWorkspace='data',
OutputWorkspace='detectors')
CORELLI#
The same garnet crystal.
### experiment data info ###
IPTS = 31429
run = 324246
### -------------------- ###
### sample info ###
a = 11.93
b = 11.93
c = 11.93
alpha = 90
beta = 90
gamma = 90
centering = 'I'
### ----------- ###
### incident wavelength ###
wavelength_band = [0.6, 2.5]
### ------------------- ###
filename = '/SNS/CORELLI/IPTS-{}/nexus/CORELLI_{}.nxs.h5'.format(IPTS, run)
LoadEventNexus(Filename=filename,
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='BL9:Mot:Sample:Axis1,0,1,0,1',
Axis1='BL9:Mot:Sample:Axis2,0,1,0,1',
Axis2='BL9:Mot:Sample:Axis3,0,1,0,1',
Average=True)
PreprocessDetectorsToMD(InputWorkspace='data',
OutputWorkspace='detectors')
SNAP#
A corundum sapphire crystal (Al₂O₃) with space group R-3c.
### experiment data info ###
IPTS = 24179
run = 51255
### -------------------- ###
### sample info ###
a = 4.75
b = 4.75
c = 12.98
alpha = 90
beta = 90
gamma = 120
centering = 'R(obv)'
### ----------- ###
### incident wavelength ###
wavelength_band = [1, 4]
### ------------------- ###
filename = '/SNS/SNAP/IPTS-{}/nexus/SNAP_{}.nxs.h5'.format(IPTS, run)
LoadEventNexus(Filename=filename,
OutputWorkspace='data')
SetGoniometer(Workspace='data',
Axis0='BL9:Mot:omega,0,1,0,1',
Average=True)
PreprocessDetectorsToMD(InputWorkspace='data',
OutputWorkspace='detectors')
Converting data to Q-sample#
After loading data and seting the goniometer, a first step could be to convert to a three-dimensional Q-sample representation, the reference frame attaced to the sample. The maximum Q-sample extent is determined by the (minimum) incident wavelength and largest reachable scattering angle:
This is the same as Q-laboratory frame.
Monochromatic beam#
For data collected at HFIR, the detector counts that varies with goniometer setting is transformed into reciprocal space with fixed wavelength.
two_theta = mtd['data'].getExperimentInfo(0).run().getProperty('TwoTheta').value
Q_max = 4*np.pi/wavelength*np.sin(0.5*max(two_theta))
ConvertHFIRSCDtoMDE(InputWorkspace='data',
Wavelength=wavelength,
LorentzCorrection=True,
MinValues=[-Q_max,-Q_max,-Q_max],
MaxValues=[+Q_max,+Q_max,+Q_max],
OutputWorkspace='md')
White beam#
The time-of-flight SNS data is converted to reciprocal space with the detector coordinates.
two_theta = mtd['detectors'].column('TwoTheta')
Q_max = 4*np.pi/min(wavelength_band)*np.sin(0.5*max(two_theta))
ConvertToMD(InputWorkspace='data',
QDimensions='Q3D',
dEAnalysisMode='Elastic',
Q3DFrames='Q_sample',
LorentzCorrection=True,
MinValues=[-Q_max,-Q_max,-Q_max],
MaxValues=[+Q_max,+Q_max,+Q_max],
PreprocDetectorsWS='detectors',
OutputWorkspace='md')
It is often advised to apply the LorentzCorrection option when initially viewing data or finding peaks. This helps reduce the sloping background at low \(|Q|\).
Note
When binning data with normalization, the Lorentz correction is done as part of normalizing procedure. This option should not be used on the 3d Q-sample event MDWorkspace in these workflows.
Determining the UB-matrix#
After converting data to \(Q\)-sample reference frame, it is possible to harvest strong single crystal peaks and using them to determine an intial \(UB\)-matrix. The \(UB\)-matrix encodes the orientation of the sample (\(U\)) and the parameters of the unit cell (\(B\)) and maps the index of the peaks to momentum transfer (scattering vector).
Finding peaks#
A first strategy to determine an initial \(UB\)-matrix from a \(Q\)-sample event MDWorkspace is to find strong peaks from the data. The algorithm FindPeaksMD can be used to find single crystal peaks using three control variables;
MaxPeaks: maximum number of peaks to find
PeakDistanceThreshold: minimum \(Q\)-spacing (Å⁻¹) between peaks
DensityThresholdFactor: minimum signal per volume to be considered a peak
The best number of peaks to find is sample-dependent but a good starting number could be 50. In addition, the mininum \(Q\)-spacing is given by the largest interplanar-spacing (typically \(2\pi/d_\mathrm{max}\)). The minimum density also depends on the experiment, but a good starting vaule could be between 100 and 1000.
FindPeaksMD(InputWorkspace='md',
MaxPeaks=max_peaks,
PeakDistanceThreshold=Q_min,
DensityThresholdFactor=strong_threshold,
OutputWorkspace='strong')
Although this algorithm is useful for finding peaks, it can find false positives that include powder lines and detector edges. It also works best if the Lorentz correction is applied first.
Determining initial Niggli cell#
One way of determining the initial \(UB\)-matrix from a PeaksWorkspace is finding the Niggli (primitive) cell. This corresponds to the reduced cell with small edges that provides a unique description independent of the lattice symmetry. It is necessary to provide a lower and upper bound of lattice constants \(a\), \(b\), and \(c\) as MinD and MaxD, respectively.
FindUBUsingFFT(PeaksWorkspace='strong',
MinD=d_min,
MaxD=d_max,
Tolerance=tol)
This algorithm performs a search for lattice vectors that index the found peaks using the fast Fourier transform.
Transforming to conventional cell#
Once a reduced cell is found, it is desirable to obtain a conventional cell for further refinement. Possible form number options corresponding to a particuar Barvais lattice can be found using the ShowPossibleCells algorithm.
MaxScalarError: limits the cell choices
with capture_logs(level='notice') as logs:
ShowPossibleCells(PeaksWorkspace='strong',
MaxScalarError=max_error,
AllowPermuations=True,
BestOnly=False)
vals = logs.getvalue()
vals = [val for val in vals.split('\n') if val.startswith('Form')]
The captured logs are placed into a list of form numbers. To apply the transform to a given form number, it is necessary to use SelectCellWithForm.
SelectCellWithForm(PeaksWorkspace='strong',
FormNumber=number,
Apply=True,
Tolerance=tol)
There are 44 form numbers corresponding to various reduced cell lattice parameters [1].
|
|
|
|
|
|
Determining conventional cell directly#
Another possible way of determining the \(UB\)-matrix from a PeaksWorkspace is directly from prior known lattice parameters serving as the initial guess.
FindUBUsingLatticeParameters(PeaksWorkspace='strong',
a=a,
b=b,
c=c,
alpha=alpha,
beta=beta,
gamma=gamma,
Tolerance=tol)
Indexing peaks#
Having determined the initial lattice parameters and orientation, the peaks can be indexed within some provided tolerance.
IndexPeaks(PeaksWorkspace='strong',
Tolerance=tol,
RoundHKLs=True)
The \(UB\)-matrix is refined as a part of this process to obtain the uncertainties of the lattice parameters. For that reason, enough peaks must be
Transforming hkls#
Considering symmetry, there can be several ambigous ways of defining the sample orientation. In addition, sometimes it is desirable to change the lattice parameters and corresponding orientation completely. For these reasons, it is often desirable to transform the oriented lattice and re-index the peaks in a new setting. Definint a transformation matrix \(T\), the \(UB\) and corresponding indexing transforms according to;
The transformation is required to be invertible. In addition, the determinant should be greater than zero to preserve a right-handed coordinate system.
TransformHKL(PeaksWorkspace='strong',
Tolerance=tol,
HKLTransform=T)
Refining the UB-matrix#
Strong reflections are useful for obtaining initial orientation and lattice parameters. More peaks are often needed to check the \(UB\)-matrix and optimize it further.
Predicting peaks#
To test the initial result and further refine (constrain) the lattice parameters, it is often useful to predict the location of peak positions given the lattice centering. The algorithm PredictPeaks uses the UB-matrix defined on the input workspaces and predicts peaks with the reflection condition corresponding to the lattice centering. It is necessary to provide a minimum and maximum \(d\)-spacing and wavelength \(\lambda\);
ReflectionCondition: lattice centering
refl_cond_dict = {'P': 'Primitive',
'I': 'Body centred',
'F': 'All-face centred',
'R': 'Primitive', # rhombohedral axes
'R(obv)': 'Rhombohderally centred, obverse', # hexagonal axes
'R(rev)': 'Rhombohderally centred, reverse', # hexagonal axes
'A': 'A-face centred',
'B': 'B-face centred',
'C': 'C-face centred'}
PredictPeaks(InputWorkspace='ws',
WavelengthMin=lamda_min,
WavelengthMax=lamda_max,
MinDSpacing=d_min,
MaxDSpacing=d_max,
ReflectionCondition=refl_cond_dict[centering],
OutputWorkspace='peaks')
Symbol |
Reflection condition |
|---|---|
P |
None |
I |
\(h+k+l=2n\) |
F |
\(h,k,l\) unmixed |
R(obv) |
\(-h+k+l=3n\) |
R(rev) |
\(h-k+l=3n\) |
A |
\(k+l=2n\) |
B |
\(l+h=2n\) |
C |
\(h+k=2n\) |
Note
With all face centering, unmixed means all \(h,k,l\) are all even or odd. For trigonal space groups, R centering with rhombohedral axes has no reflection condition restrictions. However, R centering with hexagonal axes is either in the standard obverse setting (\(-h+k+l=3n\)) or reverse setting (\(h-k+l=3n\)).
Estimating peak size#
The size of a peak in reciprocal can be estimated from the strong reflections. Given a cutoff distance in reciprocal space, the algorithm PeakIntensityVsRadius uses spherical integration to construct a workspace of the intensity of each peak with increasing radius. In addition, it provides a workspace of number of peaks that satisfy a threshold signal-to-noise ratio for a given radius.
RadiusEnd: cutoff radius
NumSteps: integration points from 0 up to cutoff radius
BackgroundInnerFactor: multiplier of peak radius for inner background shell
BackgroundInnerFactor: multiplier of peak radius for outer background shell
This can be used to extract an effective radius from the assumption of a spherical peak distribution. Comparing the integrated intensity through a cutoff radius \(r\) given the standard deviation \(\sigma\), the intensity varies according to;
Illustration of cumulative spherical integration of an isotropic Gaussian distribution.#
The number of peaks at the lowest signal-to-noise threshold varying with radius is extracted using ExtractSingleSpectrum. The parameters of this fit are estimated using EstimateFitParameters. The spherical radius is typically estimated as \(3\times\sigma\) or \(4\times\sigma\).
PeakIntensityVsRadius(InputWorkspace='md',
PeaksWorkspace='strong',
RadiusEnd=r_cut,
NumSteps=25,
BackgroundInnerFactor=1,
BackgroundOuterFactor=1.5,
OutputWorkspace='intensity',
OutputWorkspace2='signal-to-noise')
ExtractSingleSpectrum(InputWorkspace='signal-to-noise',
OutputWorkspace='number-of-peaks',
WorkspaceIndex=0)
A = mtd['number-of-peaks'].extractY().max()
constraints =';constraints=(0<A<{},0<sigma<{})'.format(2*A,r_cut/3)
function = 'name=UserFunction,Formula=I*(erf(x/(sqrt(2)*sigma))-sqrt(0.5/atan(1))*(x/sigma)*exp(-0.5*(x/sigma)^2))'
EstimateFitParameters(Function=function+constraints,
InputWorkspace='number-of-peaks',
IgnoreInvalidData=True,
OutputWorkspace='params')
I, sigma = mtd['params'].column(1)
radius = 3*sigma
Refining peak positions#
Given the predicted peaks from the current-best \(UB\)-matrix, the peak positions may be slightly off due to the limited number of peaks used to define it originally. The positions in \(Q\)-sample can be refined by finding the centroid of the distribution surrounding its current position up to a given radius.
CentroidPeaksMD(InputWorkspace='md',
PeakRadius=radius,
PeaksWorkspace='peaks',
OutputWorkspace='peaks')
Refining orientation with unconstrained cell#
A refinement of the \(UB\)-matrix using indexed \(h,k,l\) and refined \(Q\)-vector positions without any restriction on the lattice parameters is given by FindUBUsingIndexedPeaks.
FindUBUsingIndexedPeaks(PeaksWorkspace='peaks',
Tolerance=tol)
Refining orientation with constrained cell#
The algorithm OptimizeLatticeForCellType can be used to refine the \(UB\)-matrix while constraining the cell to follow certain axes with restrictions on the lattice parameters.
OptimizeLatticeForCellType(PeaksWorkspace='peaks',
CellType=cell,
Apply=True,
Tolerance=tol)
Lattice system |
Lengths |
Angles |
|---|---|---|
Cubic |
\(a=b=c\) |
\(\alpha=\beta=\gamma=90^\circ\) |
Hexagonal |
\(a=b\) |
\(\alpha=\beta=90^\circ\), \(\gamma=120^\circ\) |
Rhombohedral |
\(a=b=c\) |
\(\alpha=\beta=\gamma\) |
Tetragonal |
\(a=b\) |
\(\alpha=\beta=\gamma=90^\circ\) |
Orthorhombic |
None |
\(\alpha=\beta=\gamma=90^\circ\) |
Monoclinic |
None |
\(\alpha=\gamma=90^\circ\) |
Triclinic |
None |
None |
Note
Rhombohedral refers to the axes, not the lattice centering. Also, monoclinc assumes the standard setting with \(\beta\ne90^\circ\).
Refining orientation with fixed cell#
It is also possible to determine the orientation directly while fixing all lattice constants using CalculateUMatrix. This effectively calculates the \(U\)-matrix only.
CalculateUMatrix(PeaksWorkspace='peaks',
a=a,
b=b,
c=c,
alpha=alpha,
beta=beta,
gamma=gamma)
Integrating peaks#
Once a \(UB\)-matrix is determined and refined well to index and predict the peaks, each peak can be integrated and corrected for structure refinement. Although there are various ways of integrating peaks in reciprocal space (or detector space), the main goal is to define a region-of-interest around each peak. In reciprocal space, this is typically done by defining a sphere, cylinder, or ellipsoid around the peak. A straightforward strategy involves the latter that allows the envelope size to adapt itself to the peak.
IntegratePeaksMD(InputWorkspace='md',
PeakRadius=radius,
BackgroundInnerRadius=1*radius,
BackgroundOuterRadius=1.5*radius,
PeaksWorkspace='peaks',
OutputWorkspace='integrate',
IntegrateIfOnEdge=True,
Ellipsoid=True,
FixQAxis=False,
FixMajorAxisLength=True,
AdaptiveQMultiplier=0,
UseOnePercentBackgroundCorrection=False,
UseCentroid=True,
MaxIterations=5)
Applying absorption corrections#
Although possible to apply absorption corrections to raw event data, it is often preferred to apply the corrections after integrating peaks. This allows various different correction parameters to be used such as chemical composition and those related to sample shape. The requred information is a sample shape (parameters in centimeters) and a basic description of the chemical composition. This includes the chemical formula, Z-parameter, and unit cell volume [Å].
SetSample(InputWorkspace='integrate',
Geometry={'Shape': 'FlatPlate',
'Height': 0.1,
'Width': 0.5,
'Thick': 0.25,
'Center': [0.,0.,0.]},
Material={'ChemicalFormula': 'V', # https://docs.mantidproject.org/nightly/concepts/Materials.html#materials
'ZParameter': 2.0, # must be float
'UnitCellVolume': mtd['integrate'].sample().getOrientedLattice().volume()})
AddAbsorptionWeightedPathLengths(InputWorkspace='integrate')
For single crystal diffraction experiments, the sample is rotated with a goniometer. Therefore, the sample must also rotate in order to apply the appropriate correction at each goniometer setting that the peak was measured in. The following function can be used to apply corrections to the properly oriented sample for either a ‘sphere’, ‘cylinder’, or ‘plate’.
Shape |
Horizontal |
Vertical |
Depth |
|---|---|---|---|
sphere |
diameter |
||
cylinder |
diameter |
height |
|
plate |
length |
height |
width |
Although this function provides a way to rotate the sample, the initial sample orientation must be specified. One way is to specify crystallographic directions relative to the faces of the sample. Suppose that the sample is placed in a local coordinate system with a horizontal plane and vertical direction. For the sphere, the orientation is irrelevant since the sample shape is isotropic. For the cylinder, the axis of the cylinder is aligned along the vertical direction. Similarly, for the plate, the three axes out of the face are either along the horizontal, vertical, and depth/thickness direction. The horizontal plane contains the horizontal direction and depth/thickness direction. Only two vectors are needed to describe the sample orientation, so-called \(u\)- and \(v\)-directions.
\(u\)-vector: indices of \(h\), \(k\), and \(l\) corresponding to the direction out of a face with depth (width of the sample).
\(v\)-vector: indices of \(h\), \(k\), and \(l\) corresponding to any noncollinear direction with \(u\) that is in the horizontal plane.
def AbsorptionCorrection(ws, chemical_formula, z_parameter, u_vector, v_vector, params=[0.1,0.2,0.4], shape='plate'):
"""
Apply absorption correction to a peaks table.
ws : str
Peaks table.
u_vector : list
Out-of-depth [h, k, l] indices.
v_vector : list
Indices in horizontal plane not collinear with out-of-depth vector.
params : float or list
Parameters for shape in centimeters.
shape : str
Sample shape. Either 'plate', 'cylinder', or 'sphere'.
"""
if shape == 'sphere':
assert type(params) is float
elif shape == 'cylinder':
assert len(params) == 2
else:
assert len(params) == 3
UB = mtd[ws].sample().getOrientedLattice().getUB().copy()
u = np.dot(UB, u_vector)
v = np.dot(UB, v_vector)
u /= np.linalg.norm(u)
w = np.cross(u, v)
w /= np.linalg.norm(w)
v = np.cross(w, u)
T = np.column_stack([v, w, u])
mtd[ws].run().getGoniometer().setR(T)
gamma, beta, alpha = mtd[ws].run().getGoniometer().getEulerAngles('ZYX')
if shape == 'sphere':
shape = ' \
<sphere id="sphere"> \
<radius val="{}" /> \
<centre x="0.0" y="0.0" z="0.0" /> \
<rotate x="{}" y="{}" z="{}" /> \
</sphere> \
'.format(params/100, alpha, beta, gamma)
elif shape == 'cylinder':
shape = ' \
<cylinder id="cylinder"> \
<centre-of-bottom-base r="0.0" t="0.0" p="0.0" /> \
<axis x="0.0" y="1.0" z="0" /> \
<radius val="{}" /> \
<height val="{}" /> \
<rotate x="{}" y="{}" z="{}" /> \
</cylinder> \
'.format(params[0]/100, params[1]/100, alpha, beta, gamma)
else:
shape = ' \
<cuboid id="cuboid"> \
<width val="{}" /> \
<height val="{}" /> \
<depth val="{}" /> \
<centre x="0.0" y="0.0" z="0.0" /> \
<rotate x="{}" y="{}" z="{}" /> \
</cuboid> \
<algebra val="cuboid" /> \
'.format(params[0]/100, params[1]/100, params[2]/100, alpha, beta, gamma)
shape_dict = {'Shape': 'CSG', 'Value': shape}
Rs = []
for peak in mtd[ws]:
Rs.append(peak.getGoniometerMatrix())
runs = []
for peak in mtd[ws]:
R = peak.getGoniometerMatrix()
ind = np.isclose(Rs, R).all(axis=(1,2))
i = -1 if not np.any(ind) else ind.tolist().index(True)
run = i+1
runs.append(run)
peak.setRunNumber(run)
runs = np.unique(runs).astype(int).tolist()
mat_dict = {'ChemicalFormula': chemical_formula,
'ZParameter': float(z_parameter),
'UnitCellVolume': mtd[ws].sample().getOrientedLattice().volume()}
for i, run in enumerate(runs):
FilterPeaks(InputWorkspace=ws,
FilterVariable='RunNumber',
FilterValue=run,
Operator='=',
OutputWorkspace='tmp')
R = mtd['tmp'].getPeak(0).getGoniometerMatrix()
mtd['tmp'].run().getGoniometer().setR(R)
omega, chi, phi = mtd['tmp'].run().getGoniometer().getEulerAngles('YZY')
SetGoniometer(Workspace='tmp',
Axis0='{},0,1,0,1'.format(omega),
Axis1='{},0,0,1,1'.format(chi),
Axis2='{},0,1,0,1'.format(phi))
SetSample(InputWorkspace='tmp',
Geometry=shape_dict,
Material=mat_dict)
AddAbsorptionWeightedPathLengths(InputWorkspace='tmp', ApplyCorrection=True)
if i == 0:
CloneWorkspace(InputWorkspace='tmp', OutputWorkspace=ws+'_corr')
else:
CombinePeaksWorkspaces(LHSWorkspace=ws+'_corr', RHSWorkspace='tmp', OutputWorkspace=ws+'_corr')
hkl = np.eye(3)
Q = np.matmul(UB, hkl)
reciprocal_lattice = np.matmul(R, Q)
shape = mtd['tmp'].sample().getShape()
mesh = shape.getMesh()*100
mesh_polygon = Poly3DCollection(mesh, edgecolors='k', facecolors='w', alpha=0.5, linewidths=1)
fig, ax = plt.subplots(subplot_kw={'projection': 'mantid3d', 'proj_type': 'persp'})
ax.add_collection3d(mesh_polygon)
ax.set_title('Run #{}'.format(run))
ax.set_xlabel('x [cm]')
ax.set_ylabel('y [cm]')
ax.set_zlabel('z [cm]')
ax.set_mesh_axes_equal(mesh)
ax.set_box_aspect((1,1,1))
colors = ['r', 'g', 'b']
origin = (ax.get_xlim3d()[1], ax.get_ylim3d()[1], ax.get_zlim3d()[1])
origin = (0, 0, 0)
lims = ax.get_xlim3d()
factor = (lims[1]-lims[0])/3
for j in range(3):
vector = reciprocal_lattice[:,j]
vector_norm = vector/np.linalg.norm(vector)
ax.quiver(origin[0], origin[1], origin[2],
vector_norm[0], vector_norm[1], vector_norm[2],
length=factor,
color=colors[j],
linestyle='-')
ax.view_init(vertical_axis='y', elev=27, azim=50)
fig.show()
Saving integrated intensity files#
After the peaks are integrated, there are several options for saving the PeaksWorkspace into one of several formats. This includes,
SHELX
FullProf
Jana
GSAS-II
WinGX
LAUENORM
Direction cosines [2] can be calculated for each peak and output to an intensity file using SaveHKL and SaveHKLCW.
Normalizing data#
A \(UB\)-matrix also allows for the projection of data into a reciprocal space volume with appropriate statistical corrections [3].
Making slices and cuts from histogram data#
Data normalized and projected into reciprocal space are binned uniformly along each dimension. It is possible to extract and integrate slices/cuts from the data using IntegrateMDHistoWorkspace. For each dimension projection P\(i\) (e.g. \(i=1,2,3\)), the operation along it can be parameterized by the follwing convention;
Two values (min,max);
Sum all values between minumum and maximum.
Results in one bin.
Three values with middle zero (min,0,max);
Crop between minumum and maximum.
Bin step remains the same.
Note
Normalizing data with appropriate statistical weight may result in nan values which indicate regions without any experimental measurement statistics. Summing nan values will propagate (contaminate) the resulting sums. The statistically approprate way to handle this is to perform the integration over the raw data and normalization weights separately, then divide [4]. This would be equivalent of re-processing the normalization data with larger bins.
IntegrateMDHistoWorkspace(InputWorkspace='data',
P1Bin='-1,0,2,
P2Bin='-2,0,3,
P3Bin='-1.5,0,2,
OutputWorkspace='crop')
IntegrateMDHistoWorkspace(InputWorkspace='data',
P1Bin='-1,0,2,
P2Bin='-2,0,3,
P3Bin='-1.5,2,
OutputWorkspace='slice')
IntegrateMDHistoWorkspace(InputWorkspace='data',
P1Bin='-1,2,
P2Bin='-2,3,
P3Bin='-1.5,0,2,
OutputWorkspace='cut')
Saving histogram data to ASCII#
This function can be used and/or modified to save an MDHistoWorkspace to text-based file. The first column(s) of the output corresponds to a flattened dimension. The last two columns contain the signal and uncertainty. The output format can be modified, if needed.
def SaveMDToAscii(workspace, filename, exclude_integrated=True, format='%.6e'):
"""
Save an MDHistoToWorkspace to an ASCII file (column format).
workspace : str
Name of workspace as string.
filename : str
Path to output file.
exclude_integrated : bool, optional
Exclude integrated dimensions with bin size of one. Default is `True`.
format : str
Column format.
"""
ws = mtd[workspace]
if ws.id() != 'MDHistoWorkspace':
raise ValueError('The workspace is not an MDHistoToWorkspace')
if exclude_integrated:
dims = ws.getNonIntegratedDimensions()
else:
dims = [ws.getDimension(i) for i in range(ws.getNumDims())]
dimarrays = [np.arange(d.getMinimum() + (d.getX(1) - d.getX(0)) * 0.5, \
d.getMaximum(), d.getX(1) - d.getX(0)) for d in dims]
if len(dimarrays) > 1:
newdimarrays = np.meshgrid(*dimarrays, indexing='ij')
else:
newdimarrays = dimarrays
data = ws.getSignalArray()*1.
err = np.sqrt(ws.getErrorSquaredArray())
header = 'Intensity Error '+' '.join([d.getName() for d in dims])
header += '\n shape: '+'x'.join([str(d.getNBins()) for d in dims])
to_save = np.c_[data.flatten(), err.flatten()]
for d in newdimarrays:
to_save = np.c_[to_save, d.flatten()]
np.savetxt(filename, to_save, fmt=format, header=header)
SaveMDToAscii('data', '/path/to/output/file.txt')