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.

        %%{init: {'theme':'forest'}}%%

flowchart TB
    spice(DEMAND \n shared/autoreduce/*.nxs data) -->|<a href="https://docs.mantidproject.org/nightly/algorithms/HB3AAdjustSampleNorm-v1.html">HB3AAdjustSampleNorm</a>| histo(Detector image vs rotation index <a href="https://docs.mantidproject.org/nightly/concepts/MDHistoWorkspace.html">MDHistoWorkspace</a>)
    epics(WAND<sup>2</sup> \n nexus/*.nxs.h5 data) -->|<a href="https://docs.mantidproject.org/nightly/algorithms/LoadWANDSCD-v1.html">LoadWANDSCD</a>| histo
    sns[MANDI/TOPAZ/CORELLI/SNAP \n nexus/*.nxs.h5 data] -->|<a href="https://docs.mantidproject.org/nightly/algorithms/LoadEventNexus-v1.html">LoadEventNexus</a>| event(TOF vs counts <a href="https://docs.mantidproject.org/nightly/concepts/EventWorkspace.html">EventWorkspace</a>)
    histo -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SetGoniometer-v1.html">SetGoniometer</a> <a href="https://docs.mantidproject.org/nightly/algorithms/ConvertHFIRSCDToMDE-v1.html">ConvertHFIRSCDToMDE</a>| md[Q-sample <a href="https://docs.mantidproject.org/nightly/concepts/MDWorkspace.html">MDEventWorkspace</a>]
    event -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SetGoniometer-v1.html">SetGoniometer</a> <a href="https://docs.mantidproject.org/nightly/algorithms/ConvertToMD-v1.html">ConvertToMD</a>| md
    

WAND²

A salt crystal example (NaCl) described by space group Fm-3m.

Example loading several WAND² data as an MDHistoWorkspace.

### experiment data info ###
IPTS = 22745
runs = range(1461453, 1463252)
### -------------------- ###

### 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.

Example loading DEMAND scan data as an MDHistoWorkspace.

### 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.

Example loading MANDI data as an EventWorkspace.

### 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.

Example loading TOPAZ data as an EventWorkspace.

### 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.

Example loading CORELLI data as an EventWorkspace.

### 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.

Example loading SNAP data as an EventWorkspace.

### 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:

\[Q_\mathrm{max}=\frac{4\pi}{\lambda_\mathrm{min}}\sin{\theta_\mathrm{max}}\]

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.

Example converting MDHistoWorkspace to 3d Q-sample event MDWorkspace with ConvertHFIRSCDtoMDE.

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.

Example converting EventWorkspace to 3d Q-sample event MDWorkspace with ConvertToMD.

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).

        %%{init: {'theme':'forest'}}%%
flowchart TB
    md[Q-sample <br> <a href="https://docs.mantidproject.org/nightly/concepts/MDWorkspace.html">MDEventWorkspace</a>] 
    md --> |<a href="https://docs.mantidproject.org/nightly/algorithms/FindPeaksMD-v1.html">FindPeaksMD</a>| strong(Strong <br> <a href="https://docs.mantidproject.org/nightly/concepts/PeaksWorkspace.html">PeaksWorkspace</a>)
    strong --> |<a href="https://docs.mantidproject.org/nightly/algorithms/FindUBUsingFFT-v1.html">FindUBUsingFFT</a> <br> <a href="https://docs.mantidproject.org/nightly/algorithms/IndexPeaks-v1.html">IndexPeaks</a>| primitive[Niggli <br> <a href="https://docs.mantidproject.org/nightly/concepts/Lattice.html">Cell</a>]
    
    strong --> |<a href="https://docs.mantidproject.org/nightly/algorithms/FindUBUsingLatticeParameters-v1.html">FindUBUsingLatticeParameters</a> <br> <a href="https://docs.mantidproject.org/nightly/algorithms/IndexPeaks-v1.html">IndexPeaks</a>| conventional[Conventional <br> <a href="https://docs.mantidproject.org/nightly/concepts/Lattice.html">Cell</a>]
    
    primitive --> |<a href="https://docs.mantidproject.org/nightly/algorithms/ShowPossibleCells-v1.html">ShowPossibleCells</a> <br> <a href="https://docs.mantidproject.org/nightly/algorithms/SelectCellWithForm-v1.html">SelectCellWithForm</a>| conventional
    
    conventional --> |<a href="https://docs.mantidproject.org/nightly/algorithms/PredictPeaks-v1.html">PredictPeaks</a> <br> <a href="https://docs.mantidproject.org/nightly/algorithms/CentroidPeaksMD-v1.html">CentroidPeaksMD</a> <br> <a href="https://docs.mantidproject.org/nightly/algorithms/IntegratePeaksMD-v2.html">IntegratePeaksMD</a>| refined(Refined <br> <a href="https://docs.mantidproject.org/nightly/concepts/PeaksWorkspace.html">PeaksWorkspace</a>)
    
    refined --> |<a href="https://docs.mantidproject.org/nightly/algorithms/FindUBUsingIndexedPeaks-v1.html">FindUBUsingIndexedPeaks</a>| unconstrained[Unconstrained <br> <a href="https://docs.mantidproject.org/nightly/concepts/Lattice.html">Cell</a>]
    
    refined --> |<a href="https://docs.mantidproject.org/nightly/algorithms/OptimizeLatticeForCellType-v1.html">OptimizeLatticeForCellType</a>| constrained[Constrained <br> <a href="https://docs.mantidproject.org/nightly/concepts/Lattice.html">Cell</a>]
    
    refined --> |<a href="https://docs.mantidproject.org/nightly/algorithms/CalculateUMatrix-v1.html">CalculateUMatrix</a>| fixed[Fixed <br> <a href="https://docs.mantidproject.org/nightly/concepts/Lattice.html">Cell</a>]
    
    unconstrained --> |<a href="https://docs.mantidproject.org/nightly/algorithms/SaveIsawUB-v1.html">SaveIsawUB</a>| ub(UB .mat file)
    constrained --> |<a href="https://docs.mantidproject.org/nightly/algorithms/SaveIsawUB-v1.html">SaveIsawUB</a>| ub
    fixed --> |<a href="https://docs.mantidproject.org/nightly/algorithms/SaveIsawUB-v1.html">SaveIsawUB</a>| ub


    

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.

Finding strong peaks with FindPeaksMD.

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.

• MinD: lower limit of \(a\), \(b\), and \(c\)

• MaxD: upper limit of \(a\), \(b\), and \(c\)

Determining oriented Niggli cell with FindUBUsingFFT.

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 based on the calculated errors

Listing possible conventional cells with ShowPossibleCells.

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.

• FormNumber: 1-44 redued form numbers

Choosing the convetional cell with SelectCellWithForm.

SelectCellWithForm(PeaksWorkspace='strong',
                   FormNumber=number,
                   Apply=True,
                   Tolerance=tol)

There are 44 form numbers corresponding to various reduced cell lattice parameters [1].

Reduced form numbers.

\(a=b=c\) #1 Cubic F #2 Rhombohedral R #3 Cubic P #4 Rhombohedral R #5 Cubic I #6 Tetragonal I #7 Tetragonal I #8 Orthorhombic I

\(a=b\) #9 Rhombohedral R #10 Monoclinic C/I #11 Tetragonal P #12 Hexagonal P #13 Orthorhombic C #14 Monoclinic C/I #15 Tetragonal I #16 Orthorhombic F #17 Monoclinic I/C

\(b=c\) #18 Tetragonal I #19 Orthorhombic I #20 Monoclinic C/I #21 Tetragonal P #22 Hexagonal P #23 Orthorhombic C #24 Rhombohedral R #25 Monoclinic C/I

\(a\le b\le c\)

#26 Orthorhombic F #27 Monoclinic I/C #28 Monoclinic C #29 Monoclinic C #30 Monoclinic C #31 Triclinic P #32 Orthorhombic P

#33 Monoclinic P #34 Monoclinic P #35 Monoclinic P #36 Orthorhombic C #37 Monoclinic C/I #38 Orthorhombic C

#39 Monoclinic C/I #40 Orthorhombic C #41 Monoclinic C/I #42 Orthorhombic I #43 Monoclinic I #44 Triclinic P

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.

• a, b, c: \(a\), \(b\), \(c\)-lattice constants (angstroms)

• alpha, beta, gamma: \(\alpha\), \(\beta\), \(\gamma\)-lattice constants (degrees)

Determining oriented conventional cell with FindUBUsingLatticeParameters.

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.

Indexing peaks with found \(UB\)-matrix using IndexPeaks.

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;

\[\begin{split}\begin{align} \begin{bmatrix} UB_{11}^\prime & UB_{12}^\prime & UB_{13}^\prime \\ UB_{21}^\prime & UB_{22}^\prime & UB_{23}^\prime \\ UB_{31}^\prime & UB_{32}^\prime & UB_{33}^\prime \\ \end{bmatrix}&= \begin{bmatrix} T_{11} & T_{12} & T_{13} \\ T_{21} & T_{22} & T_{23} \\ T_{31} & T_{32} & T_{33} \\ \end{bmatrix}^{-1} \begin{bmatrix} UB_{11} & UB_{12} & UB_{13} \\ UB_{21} & UB_{22} & UB_{23} \\ UB_{31} & UB_{32} & UB_{33} \\ \end{bmatrix} \\ \begin{bmatrix} h^\prime \\ k^\prime \\ l^\prime \end{bmatrix}&= \begin{bmatrix} T_{11} & T_{12} & T_{13} \\ T_{21} & T_{22} & T_{23} \\ T_{31} & T_{32} & T_{33} \\ \end{bmatrix} \begin{bmatrix} h \\ k \\ l \end{bmatrix}\\ \end{align}\end{split}\]

The transformation is required to be invertible. In addition, the determinant should be greater than zero to preserve a right-handed coordinate system.

Transforming peaks and oriented lattice with TransformHKL.

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

Predicting peaks with PredictPeaks.

refl_cond_dict = {'P': 'Primitive',
                  'I': 'Body centred',
                  'F': 'All-face centred',
                  'R': 'Primitive', # rhombohedral axes
                  'R(obv)': 'Rhombohedrally centred, obverse', # hexagonal axes
                  'R(rev)': 'Rhombohedrally 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')

Integral reflection conditions for centered cells.

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, with rhombohedral axes has no reflection condition lattice centering restrictions similar to a “primitive” lattice. 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;

\[I \bigg[ \mathrm{erf}\,\bigg(\frac{r}{\sqrt{2}\sigma}\bigg)-\sqrt{\frac{2}{\pi}}\frac{r}{\sigma}\exp\bigg(-\frac{r^2}{2\sigma^2}\bigg) \bigg]\]

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\).

Estimating peak size with PeakIntensityVsRadius.

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=A*(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.

Refining peak locations with CentroidPeaksMD.

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.

Refining UB with 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.

Refining UB with OptimizeLatticeForCellType.

OptimizeLatticeForCellType(PeaksWorkspace='peaks',
                           CellType=cell,
                           Apply=True,
                           Tolerance=tol)

Parameter constraints for each lattice system.

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.

Refining UB with CalculateUMatrix.

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.

Integrating peak intensities with IntegratePeaksMD.

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)

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

        %%{init: {'theme':'forest'}}%%

flowchart TB
    peaks[Integrated <a href="https://docs.mantidproject.org/nightly/concepts/PeaksWorkspace.html">PeaksWorkspace</a>] -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SaveHKL-v1.html">SaveHKL</a>| hkl(Shelx/GSAS/Jana *.hkl)
    peaks -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SaveReflections-v1.html">SaveReflections</a>| int(FullProf *.int)
    peaks -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SaveHKLCW-v1.html">SaveHKLCW</a>| wingx(WinGX *.hkl)
    peaks -->|<a href="https://docs.mantidproject.org/nightly/algorithms/SaveLauenorm-v1.html">SaveLauenorm</a>| mtz(LAUENORM laue)
    

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].

        %%{init: {'theme':'forest'}}%%

flowchart TB
    vanmat[Vanadium \n *.nxs data] --> |<a href="https://docs.mantidproject.org/nightly/algorithms/LoadNexus-v1.html">LoadNexus</a>| spectrum(Solid angle and flux \n <a href="https://docs.mantidproject.org/nightly/concepts/MatrixWorkspace.html">MatrixWorkspace</a>)
    vanhist[Vanadium \n *.nxs data] --> |<a href="https://docs.mantidproject.org/nightly/algorithms/LoadMD-v1.html">LoadMD</a>| counts(Detector counts \n <a href="https://docs.mantidproject.org/nightly/concepts/MDHistoWorkspace.html">MDHistoWorkspace</a>)
    epics(WAND<sup>2</sup> \n nexus/*.nxs.h5 data) --> |<a href="https://docs.mantidproject.org/nightly/algorithms/LoadWANDSCD-v1.html">LoadWANDSCD</a>| histo
    sns[MANDI/TOPAZ/CORELLI/SNAP \n nexus/*.nxs.h5 data] --> |<a href="https://docs.mantidproject.org/nightly/algorithms/LoadEventNexus-v1.html">LoadEventNexus</a>| event(TOF vs counts \n <a href="https://docs.mantidproject.org/nightly/concepts/EventWorkspace.html">EventWorkspace</a>)
    event --> |<a href="https://docs.mantidproject.org/nightly/algorithms/SetGoniometer-v1.html">SetGoniometer</a> \n <a href="https://docs.mantidproject.org/nightly/algorithms/ConvertToMD-v1.html">ConvertToMD</a>| mde[Q-sample \n <a href="https://docs.mantidproject.org/nightly/concepts/MDWorkspace.html">MDEventWorkspace</a>]
    spice(DEMAND \n shared/autoreduce/*.nxs data) -->|<a href="https://docs.mantidproject.org/nightly/algorithms/HB3AAdjustSampleNorm-v1.html">HB3AAdjustSampleNorm</a>| histo(Detector image vs rotation index \n <a href="https://docs.mantidproject.org/nightly/concepts/MDHistoWorkspace.html">MDHistoWorkspace</a>)
    mde --- mdnorm[ ]:::empty
    spectrum --- mdnorm
    mdnorm --> |<a href="https://docs.mantidproject.org/nightly/algorithms/MDNorm-v1.html">MDNorm</a>| mdh[HKL/Q-sample \n <a href="https://docs.mantidproject.org/nightly/concepts/MDWorkspace.html">MDHistoWorkspace</a>]
    histo --- |<a href="https://docs.mantidproject.org/nightly/algorithms/SetGoniometer-v1.html">SetGoniometer</a>| convertq[ ]:::empty
    counts --- convertq
    convertq --> |<a href="https://docs.mantidproject.org/nightly/algorithms/ConvertWANDSCDtoQ-v1.html">ConvertWANDSCDtoQ</a>| mdh
    classDef empty width:-1px,height:-1px;
    

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;

1. Two values (min,max);

   • Sum all values between minumum and maximum.

   • Results in one bin.

2. 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.

Example croping a volume into a smaller volume. Bin size remains the same in all dimensions.

IntegrateMDHistoWorkspace(InputWorkspace='data',
                          P1Bin='-1,0,2,
                          P2Bin='-2,0,3,
                          P3Bin='-1.5,0,2,
                          OutputWorkspace='crop')

Example slicing a volume into a two-dimensional planar slice. The third dimension is integrated.

IntegrateMDHistoWorkspace(InputWorkspace='data',
                          P1Bin='-1,0,2,
                          P2Bin='-2,0,3,
                          P3Bin='-1.5,2,
                          OutputWorkspace='slice')

Example cutting a volume into a one-dimensional linear cut. The first and second dimensions are integrated.

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.

Example saving MDHistoWorkspace to ASCII file.

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')

References

[1]

Mighell, A. D. Lattice Symmetry and Identification-The Fundamental Role of Reduced Cells in Materials Characterization. J Res Natl Inst Stand Technol 2001, 106 (6), 983–995. https://doi.org/10.6028/jres.106.050.

[2]

Katrusiak, A. Absorption Correction for Crystal-Environment Attachments from Direction Cosines. Zeitschrift für Kristallographie - Crystalline Materials 2001, 216 (12), 646–647. https://doi.org/10.1524/zkri.216.12.646.22488.

[3]

Michels-Clark, T. M.; Savici, A. T.; Lynch, V. E.; Wang, X.; Hoffmann, C. M. Expanding Lorentz and Spectrum Corrections to Large Volumes of Reciprocal Space for Single-Crystal Time-of-Flight Neutron Diffraction. J Appl Crystallogr 2016, 49 (Pt 2), 497–506. https://doi.org/10.1107/S1600576716001369.

[4]

Savici, A. T.; Gigg, M. A.; Arnold, O.; Tolchenov, R.; Whitfield, R. E.; Hahn, S. E.; Zhou, W.; Zaliznyak, I. A. Efficient Data Reduction for Time-of-Flight Neutron Scattering Experiments on Single Crystals. J Appl Cryst 2022, 55 (6), 1514–1527. https://doi.org/10.1107/S1600576722009645.