diff --git a/documentation/user_docs/docs/manual/Correcting_for_sample_misalignment.rst b/documentation/user_docs/docs/manual/Correcting_for_sample_misalignment.rst index 6358f5b819..8c229b46f5 100644 --- a/documentation/user_docs/docs/manual/Correcting_for_sample_misalignment.rst +++ b/documentation/user_docs/docs/manual/Correcting_for_sample_misalignment.rst @@ -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: @@ -35,7 +36,7 @@ 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, ) @@ -43,9 +44,9 @@ from the expected values. 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}\}`. @@ -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 @@ -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:: @@ -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, ); + alignment_info = refine_crystal(rlu_actual, alatt, angdeg, rlu_expected, ); 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: @@ -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 ----------------------------------------- @@ -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): @@ -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:: @@ -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) ======================================================================= @@ -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, ) + [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, ) The inputs are: @@ -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: diff --git a/horace_core/lattice_functions/bragg_positions.m b/horace_core/lattice_functions/bragg_positions.m index fb7e9a76f0..9e10ae00c9 100644 --- a/horace_core/lattice_functions/bragg_positions.m +++ b/horace_core/lattice_functions/bragg_positions.m @@ -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 @@ -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] % @@ -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 @@ -111,7 +112,7 @@ 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{:}); @@ -119,8 +120,8 @@ % 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); @@ -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 @@ -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 diff --git a/horace_core/lattice_functions/refine_crystal.m b/horace_core/lattice_functions/refine_crystal.m index 09ab70d10a..32e1086e81 100644 --- a/horace_core/lattice_functions/refine_crystal.m +++ b/horace_core/lattice_functions/refine_crystal.m @@ -1,8 +1,8 @@ -function [alignment_info, error_estimate] = refine_crystal(rlu0,alatt0,angdeg0,rlu,varargin) +function [alignment_info, error_estimate] = refine_crystal(rlu_actual,alatt0,angdeg0,rlu_expected,varargin) % Refine crystal orientation and lattice parameters % -% >> alignment_info = refine_crystal(rlu0, alatt0, angdeg0, rlu) -% >> alignment_info = refine_crystal(rlu0, alatt0, angdeg0, rlu, alatt_init, angdeg_init) +% >> alignment_info = refine_crystal(rlu_actual, alatt0, angdeg0, rlu_expected) +% >> alignment_info = refine_crystal(rlu_actual, alatt0, angdeg0, rlu_expected, alatt_init, angdeg_init) % Return estimates of the error in the crystal alignment fit % >> [..., error_estimate] = refine_crystal( ... ) % @@ -22,13 +22,14 @@ % % Input: % ------ -% rlu0 Positions of reciprocal lattice vectors as h,k,l in reference lattice +% rlu_actual Positions of reciprocal lattice vectors as h,k,l in reference lattice % (n x 3 matrix, n=no. reflections) % alatt0 Reference lattice parameters [a,b,c] (Angstroms), which -% would provide bragg peaks at rlu0 +% would provide bragg peaks at rlu_actual % angdeg0 Reference lattice angles [alph,bet,gam] (deg), which -% would provide bragg peaks at rlu0 -% rlu True indexes of reciprocal lattice vectors (n x 3 matrix) +% would provide bragg peaks at rlu_actual +% rlu_expected True indices (the indices would be expected in aligned lattice) +% of reciprocal lattice vectors (n x 3 matrix) % % Optional input parameter: % alatt_init Initial lattice parameters for start of refinement [a,b,c] (Angstroms) @@ -74,7 +75,7 @@ % in the two frames are related by % v(i)= rotmat(i,j)*v0(j) % distance Distances between peak positions and points given by true indexes, in input -% argument rlu, in the refined crystal lattice. (Ang^-1) +% argument rlu_expected, in the refined crystal lattice. (Ang^-1) % % rotangle Angle of rotation corresponding to rotmat (to give a measure % of the misorientation) (degrees) @@ -87,11 +88,11 @@ % % EXAMPLES % Want to refine crystal orientation only: -% >> alignment_info =refine_crystal (rlu0, alatt0, angdeg0, rlu, 'fix_lattice') +% >> alignment_info =refine_crystal (rlu_actual, alatt0, angdeg0, rlu_expected, 'fix_lattice') % the alignment info would contain the initial lattice values i.e. alatt0, angdeg0 % % Want to refine lattice parameters a,b,c as well as crystal orientation: -% >> alignment_info=refine_crystal (rlu0, alatt0, angdeg0, rlu, 'fix_angdeg') +% >> alignment_info=refine_crystal (rlu_actual, alatt0, angdeg0, rlu_expected, 'fix_angdeg') % the alignment info would contain the initial lattice angles i.e. angdeg0 small=1e-10; @@ -105,7 +106,7 @@ check_options_consistency(present,opt); % Check input arguments -lattice0 = check_input_arguments(rlu0,rlu,alatt0,angdeg0); +lattice0 = check_input_arguments(rlu_actual,rlu_expected,alatt0,angdeg0); % Check initial lattice parameters for refinement, if given, are acceptable lattice_init = check_additional_args(lattice0,args{:}); @@ -116,8 +117,8 @@ b0 = bmatrix(lattice0(1:3),lattice0(4:6)); binit = bmatrix(lattice_init(1:3),lattice_init(4:6)); -vcryst0=b0*rlu0'; % crystal Cartesian coords in reference lattice -vcryst_init=binit*rlu'; % crystal Cartesian coords in initial lattice +vcryst0=b0*rlu_actual'; % crystal Cartesian coords in reference lattice +vcryst_init=binit*rlu_expected'; % crystal Cartesian coords in initial lattice % Check lengths are all non-zero lensqr0=sum(vcryst0.^2,1); @@ -136,7 +137,7 @@ % Fit % --- -nv=size(rlu,1); +nv=size(rlu_expected,1); vcryst0_3=repmat(vcryst0',3,1); % Treble the number of vectors, as will compute three components of deviations pars=[lattice_init,rotvec_ave(:)']; @@ -166,7 +167,7 @@ end kk = multifit (vcryst0_3, zeros(3*nv,1), 0.01*ones(3*nv,1)); -kk = kk.set_fun (@reciprocal_space_deviation, {pars,rlu}, pfree, pbind); +kk = kk.set_fun (@reciprocal_space_deviation, {pars,rlu_expected}, pfree, pbind); kk = kk.set_options ('list', 0); kk = kk.set_options ('fit',[1e-4,50,-1e-6]); [distance,fitpar] = kk.fit(); @@ -192,8 +193,8 @@ distance=sqrt(sum(reshape(distance,3,nv).^2,1))'; error_estimate = struct('alatt', fitpar.sig(1:3), ... - 'angdeg', fitpar.sig(4:6), ... - 'rotvec', fitpar.sig(7:9)); + 'angdeg', fitpar.sig(4:6), ... + 'rotvec', fitpar.sig(7:9)); alignment_info = crystal_alignment_info(alatt,angdeg,rotvec,distance); @@ -201,10 +202,10 @@ %============================================================================================================ % Distance function for fitting %============================================================================================================ -function dist = reciprocal_space_deviation (x1,x2,x3,p,rlu) +function dist = reciprocal_space_deviation (x1,x2,x3,p,rlu_expected) % Function to calculate the distance between a point in reciprocal space and corresponding point in a reference orthonormal frame % -% >> dist = reciprocal_space_deviation (v0,p,rlu) +% >> dist = reciprocal_space_deviation (v0,p,rlu_expected) % % Input: % ------- @@ -216,14 +217,14 @@ % theta1,theta2,theta3 components of rotation vector linking % crystal Cartesian coordinates % v(i)=R_theta(i,j)*v0(j) -% rlu Components along a*, b*, c* in lattice defined by p (n x 3 array) +% rlu_expected Components along a*, b*, c* in lattice defined by p (n x 3 array) % % Output: % ------- % dist Column vector of deviations along x,y,z axes of reference crystal -% Cartesian coordinates for each of the vectors rlu in turn +% Cartesian coordinates for each of the vectors rlu_expected in turn -nv=size(rlu,1); +nv=size(rlu_expected,1); alatt=p(1:3); angdeg=p(4:6); @@ -232,7 +233,7 @@ b=bmatrix(alatt,angdeg); R=rotvec_to_rotmat2(rotvec); rlu_to_cryst0=R\b; -v=(rlu_to_cryst0*rlu')'; +v=(rlu_to_cryst0*rlu_expected')'; dv=v-[x1(1:nv),x2(1:nv),x3(1:nv)]; dist=reshape(dv',3*nv,1)./repmat(sqrt(sum(v.^2,2)),3,1); @@ -386,18 +387,18 @@ function check_options_consistency(present,opt) 'Check consistency of options to fix lattice parameters') end % -function lattice0 = check_input_arguments(rlu0,rlu,alatt0,angdeg0) -if size(rlu0,2)~=3 || size(rlu0,1)<2 || numel(size(rlu0))~=2 +function lattice0 = check_input_arguments(rlu_actual,rlu_expected,alatt0,angdeg0) +if size(rlu_actual,2)~=3 || size(rlu_actual,1)<2 || numel(size(rlu_actual))~=2 error('HORACE:lattice_functions:invalid_argument', ... 'Must be at least two input reciprocal lattice vectors, each given as triples (h,k,l)') end -if numel(size(rlu))~=2 || ~all(size(rlu0)==size(rlu)) +if numel(size(rlu_expected))~=2 || ~all(size(rlu_actual)==size(rlu_expected)) error('HORACE:lattice_functions:invalid_argument', ... 'Must be the same number of reciprocal lattice vectors in reference and new coordinate frames, each given as a triple (h,k,l)') end -if ~all(isfinite(rlu0(:))) % catch case of rlu not being found - a common input is from bragg_positions +if ~all(isfinite(rlu_actual(:))) % catch case of rlu_expected not being found - a common input is from bragg_positions error('HORACE:lattice_functions:invalid_argument', ... - 'One or more positions of the true Bragg peak positions (input argument ''rlu'') is not finite.') + 'One or more positions of the true Bragg peak positions (input argument ''rlu_expected'') is not finite.') end if isnumeric(alatt0) && numel(alatt0)==3 && all(alatt0>0) && isnumeric(angdeg0) && numel(angdeg0)==3 && all(angdeg0>0)