Skip to content

Commit

Permalink
Merge pull request #1687 from pace-neutrons/1671_alignment_doc
Browse files Browse the repository at this point in the history
alignment documentation
  • Loading branch information
abuts authored May 29, 2024
2 parents adfa845 + d09f198 commit db8435f
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 64 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,9 @@ Correcting for sample misalignment


When mounting your sample on a spectrometer, it can often be the case that it is slightly misaligned with respect to the
'perfect' alignment assumed when generating the SQW file (the ``u`` and ``v`` vectors provided in ``gen_sqw`` and
``accumulate_sqw``). It is straightforward to correct this misalignment, once enough data have been accumulated, by
'perfect' alignment assumed when generating the SQW file (the direction of ``u`` and ``v`` vectors provided in ``gen_sqw`` and
``accumulate_sqw``, where ``u`` is parallel to the beam and ``v`` defines the sample rotation plain).
It is straightforward to correct this misalignment, once enough data have been accumulated, by
comparing the positions of Bragg peaks with what they are expected to be.

Alignment correction is a two-step process:
Expand All @@ -35,17 +36,17 @@ from the expected values.

::

[rlu0, widths, wcut, wpeak] = bragg_positions (sqw, bragg_positions, ...
[rlu_actual, widths, wcut, wpeak] = bragg_positions (sqw_obj, bragg_expected, ...
radial_cut_length, radial_bin_width, radial_thickness,...
trans_cut_length, trans_bin_width, trans_thickness, ...
energy_window, <keyword options>)


The inputs are:

- ``sqw`` - the uncorrected data
- ``sqw_obj`` - ``sqw`` object with misaligned data

- ``bragg_positions`` - an n-by-3 array specifying the expected Bragg positions
- ``bragg_expected`` - an n-by-3 array specifying the Bragg positions expected from aligned crystal.

- ``radial_cut_length`` - lengths of the various cuts along each of :math:`\{\vec{Q}\}`.

Expand Down Expand Up @@ -93,7 +94,7 @@ For fitting:

The outputs are:

- ``rlu0`` - the actual peak positions as an n-by-3 matrix in :math:`h,k,l` as indexed with respect to the current
- ``rlu_actual`` - the actual peak positions as an n-by-3 matrix in :math:`h,k,l` as indexed with respect to the current
lattice parameters.

- ``widths`` - an n-by-3 array containing the FWHH in Ang^-1 of the peaks along each of the three projection axes
Expand Down Expand Up @@ -128,7 +129,7 @@ everything, using outputs from ``bragg_positions`` described above.
bragg_positions_view(wcut,wpeak)


You will be prompted in the Matlab command window as to which plot and fit you wish to view.
You will be prompted in the MATLAB command window as to which plot and fit you wish to view.

.. note::

Expand All @@ -147,21 +148,21 @@ You will be prompted in the Matlab command window as to which plot and fit you w
Step 3 - calculate the misalignment correction
----------------------------------------------

Using the outputs of ``bragg_positions``, you can determine a transformation matrix to go from the original
Using the outputs of ``bragg_positions``, you can determine a transformation to go from the original
misaligned frame to the aligned frame.

::

al_info = refine_crystal(rlu0, alatt, angdeg, bragg_peaks, <keyword options>);
alignment_info = refine_crystal(rlu_actual, alatt, angdeg, rlu_expected, <keyword options>);


The inputs are:

- ``rlu0`` - the an n-by-3 matrix of actual peak positions as in :math:`h,k,l` as indexed with the current lattice parameters
- ``rlu_actual`` - the an n-by-3 matrix of actual peak positions as in :math:`h,k,l` as indexed with the current lattice parameters

- ``alatt, angdeg`` - the lattice parameters and angles used in the original sqw file.
- ``alatt, angdeg`` - the lattice parameters and angles used in the original (misaligned) sqw file.

- ``bragg_peaks`` - the predicted (integer) Bragg peaks corresponding to ``rlu0``
- ``rlu_expected`` - the predicted (integer) Bragg peaks corresponding to ``bragg_expected``

The keyword options are:

Expand Down Expand Up @@ -190,7 +191,9 @@ The keyword options are:

To achieve finer control of the refinement of the lattice parameters, use ``free_alatt`` and ``free_angdeg``

The output is an ``crystal_alignment_info`` object which contains all the relevant data for crystal realignment.
The output is an ``crystal_alignment_info`` object which contains all the relevant data for crystal realignment, namely
the rotation matrix which aligns Crystal Cartesian frame into correct position and modified lattice parameters, if
``refine_crystal`` modified them.

Step 4 - apply the correction to the data
-----------------------------------------
Expand All @@ -203,7 +206,8 @@ There are different to do this, for different circumstances:

- When you have a loaded ``sqw`` object:

Apply the correction to the object
Apply the correction to the object. The object's lattice and pixels orientation will be modified
to aligned values.

- When you are still accumulating data (e.g. on the beamline):

Expand All @@ -219,15 +223,16 @@ There is a simple routine to apply the changes to an existing file, without the

change_crystal(win, alignment_info)

where ``alignment_info`` was determined in the steps described above. From this point out the alignment will be applied whenever pixels are loaded or manipulated (e.g. loading, cutting, plotting, etc.).
where ``alignment_info`` was determined in the steps described above. This procedure modifies lattice parameters and
the pixels will be aligned whenever they are loaded or manipulated (e.g. accessing pixel data, cutting, plotting, etc.).

Once you have confirmed that the alignment you have is the correct one, it is possible to fix the alignment to avoid this calculation step.

This is done through the ``apply_alignment`` function:
This is done through the ``finalize_alignment`` function:

::

[wout, rev_corr] = apply_alignment(win, ['-keep_original'])
[wout, rev_corr] = finalize_alignment(win, ['-keep_original'])

.. warning::

Expand All @@ -241,12 +246,16 @@ Where:

.. note::

If you use ``'-keep_original'`` you may wish to ``save`` your object as the temporary file will be cleared when the object is. (see: file_backed_objects)
If you use ``'-keep_original'`` you may wish to ``save`` your object as the temporary file will be cleared when the ``wout`` object is. (see: file_backed_objects)

- ``wout`` - Resulting ``sqw`` object or the filename to which the alignment was applied.
- ``wout`` - Resulting ``sqw`` object to which the alignment was applied. If input was kept in file or was filebacked, the object will be filebacked.

- ``rev_corr`` - A corresponding ``crystal_alignment_info`` to be able to reverse the application.

.. note::

Finalize alignment of large ``sqw`` object may take substantial time. The time may be even bigger than regenerating this file from scratch as parallel
generation is currently possible for ``sqw`` files generation but not yet implemented for ``finalize_alignment``.

Option 2 : calculate goniometer offsets for regeneration of sqw file(s)
=======================================================================
Expand All @@ -255,7 +264,7 @@ In this case there is a single routine to calculate the new goniometer offsets,

::

[alatt, angdeg, dpsi_deg, gl_deg, gs_deg] = crystal_pars_correct(u, v, alatt0, angdeg0, omega0_deg, dpsi0_deg, gl0_deg, gs0_deg, al_info, <keyword options>)
[alatt, angdeg, dpsi_deg, gl_deg, gs_deg] = crystal_pars_correct(u, v, alatt0, angdeg0, omega0_deg, dpsi0_deg, gl0_deg, gs0_deg, alignment_info, <keyword options>)


The inputs are:
Expand All @@ -277,7 +286,7 @@ The inputs are:
:math:`\text{d}\psi`, :math:`g_l` and :math:`g_s` refer to the Euler angles relative to the scattering plane. Naming
conventions may differ in other notations, e.g. :math:`\theta, \phi, \chi`.

- ``al_info`` - The ``crystal_alignment_info`` object determined above.
- ``alignment_info`` - The ``crystal_alignment_info`` object determined above.

The keywords options are:

Expand Down
31 changes: 16 additions & 15 deletions horace_core/lattice_functions/bragg_positions.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [rlu0,width,wcut,wpeak]=bragg_positions(w, rlu,...
function [rlu_actual,width,wcut,wpeak]=bragg_positions(w, rlu_expected,...
radial_cut_length, radial_bin_width, radial_thickness,...
trans_cut_length, trans_bin_width, trans_thickness, varargin)
% Get actual Bragg peak positions, given initial estimates of their positions
Expand All @@ -12,11 +12,12 @@
% Input:
% ------
% w Data source (sqw file name or sqw object)
% rlu Set of Bragg peak indices in (n x 3) matrix:
% rlu_expected Set of Bragg peak indices in (n x 3) matrix:
% h1, k1, l1
% h2, k2, l2
% : : :
% These are the indices of the Bragg peaks whose
% These are the expected (If crystal is perfectly aligned) indices
% of the Bragg peaks whose
% actual positions will be found by this function
% e.g. [1,0,0; 0,1,0; 1,1,0]
%
Expand Down Expand Up @@ -66,9 +67,9 @@
%
% Output:
% -------
% rlu0 The actual peak positions as (n x 3) matrix of h,k,l as
% rlu_actual The actual peak positions as (n x 3) matrix of h,k,l as
% indexed with the current lattice parameters
% widths Array (size (n x 3)) containing the FWHH in Ang^-1 of the
% widths Array (size (n x 3)) containing the FWHH in Ang^-1 of the
% peaks along each of the three projection axes
% wcut Array of cuts, size (n x 3), along three orthogonal
% directions through each Bragg point from which the peak
Expand Down Expand Up @@ -111,16 +112,16 @@
error('Object must be sqw type')
end

[opt,eint,absolute_binning,gau] = check_and_parse_inputs( rlu,...
[opt,eint,absolute_binning,gau] = check_and_parse_inputs( rlu_expected,...
radial_cut_length, radial_bin_width, radial_thickness,...
trans_cut_length, trans_bin_width, trans_thickness, varargin{:});


% Fit Peaks
% ---------
% Initialise output arguments
szrlu = size(rlu);
rlu0 = zeros(szrlu );
szrlu = size(rlu_expected);
rlu_actual = zeros(szrlu );
width = zeros(szrlu);
npeaks = szrlu(1);
wcut=repmat(IX_dataset_1d,npeaks,3);
Expand All @@ -133,11 +134,11 @@
% Get the matrix to convert rlu to crystal Cartesian coordinates
B = bmatrix (img.alatt, img.angdeg);

peak_problem=false(size(rlu));
peak_problem=false(size(rlu_expected));

for i=1:size(rlu,1)
for i=1:size(rlu_expected,1)
% Extract Q point through which to get three orthogonal cuts
Qrlu = rlu(i,:);
Qrlu = rlu_expected(i,:);
modQ=norm(B*Qrlu(:)); % length of Q vector in Ang^-1

% Create proj for taking three orthogonal cuts
Expand Down Expand Up @@ -214,23 +215,23 @@

% Convert peak position into r.l.u.
if all(isfinite(upos0))
rlu0(i,:)= proj.transform_img_to_hkl(upos0(:));
rlu_actual(i,:)= proj.transform_img_to_hkl(upos0(:));
else
peak_problem(i,:)=~isfinite(upos0);
rlu0(i,:)=NaN;
rlu_actual(i,:)=NaN;
end
end

disp('--------------------------------------------------------------------------------')
if any(peak_problem(:))
disp('Problems determining peak position for:')
for i=1:size(rlu,1)
for i=1:size(rlu_expected,1)
if any(peak_problem(i,:))
disp(['Peak ',num2str(i),': [',num2str(Qrlu),']',' scan(s): ',num2str(find(peak_problem(i,:)))])
end
end
disp(' ')
disp(['Total number of peaks = ',num2str(size(rlu,1))])
disp(['Total number of peaks = ',num2str(size(rlu_expected,1))])
disp('--------------------------------------------------------------------------------')
end

Expand Down
Loading

0 comments on commit db8435f

Please sign in to comment.