diff --git a/hexrd/core/fitting/calibration/laue.py b/hexrd/core/fitting/calibration/laue.py index 59f65ea7b..11030dcba 100644 --- a/hexrd/core/fitting/calibration/laue.py +++ b/hexrd/core/fitting/calibration/laue.py @@ -13,6 +13,8 @@ from hexrd.core.instrument import switch_xray_source from hexrd.core.transforms import xfcapi +from hexrd.laue.simulation import simulate_laue_pattern_on_instrument, simulate_laue_pattern_on_panel + from .abstract_grain import AbstractGrainCalibrator from .lmfit_param_handling import DEFAULT_EULER_CONVENTION @@ -174,7 +176,8 @@ def _autopick_points( # run simulation # ???: could we get this from overlays? - laue_sim = self.instr.simulate_laue_pattern( + laue_sim = simulate_laue_pattern_on_instrument( + self.instr, self.plane_data, minEnergy=self.energy_cutoffs[0], maxEnergy=self.energy_cutoffs[1], @@ -507,7 +510,8 @@ def sxcal_obj_func( for det_key, panel in instr.detectors.items(): # Simulate Laue pattern: # returns xy_det, hkls_in, angles, dspacing, energy - sim_results = panel.simulate_laue_pattern( + sim_results = simulate_laue_pattern_on_panel( + panel, [hkls_idx[det_key], bmat], minEnergy=energy_cutoffs[0], maxEnergy=energy_cutoffs[1], diff --git a/hexrd/core/instrument/detector.py b/hexrd/core/instrument/detector.py index 84df5280c..78df625cb 100644 --- a/hexrd/core/instrument/detector.py +++ b/hexrd/core/instrument/detector.py @@ -1718,161 +1718,6 @@ def simulate_rotation_series( ang_pixel_size.append(self.angularPixelSize(xys_p)) return valid_ids, valid_hkls, valid_angs, valid_xys, ang_pixel_size - def simulate_laue_pattern( - self, - crystal_data, - minEnergy=5.0, - maxEnergy=35.0, - rmat_s=None, - tvec_s=None, - grain_params=None, - beam_vec=None, - ): - """ """ - if isinstance(crystal_data, PlaneData): - plane_data = crystal_data - - # grab the expanded list of hkls from plane_data - hkls = np.hstack(plane_data.getSymHKLs()) - - # and the unit plane normals (G-vectors) in CRYSTAL FRAME - gvec_c = np.dot(plane_data.latVecOps['B'], hkls) - - # Filter out g-vectors going in the wrong direction. `gvec_to_xy()` used - # to do this, but not anymore. - to_keep = np.dot(gvec_c.T, self.bvec) <= 0 - - hkls = hkls[:, to_keep] - gvec_c = gvec_c[:, to_keep] - elif len(crystal_data) == 2: - # !!! should clean this up - hkls = np.array(crystal_data[0]) - bmat = crystal_data[1] - gvec_c = np.dot(bmat, hkls) - else: - raise RuntimeError(f'argument list not understood: {crystal_data=}') - nhkls_tot = hkls.shape[1] - - # parse energy ranges - # TODO: allow for spectrum parsing - multipleEnergyRanges = False - - if hasattr(maxEnergy, '__len__'): - if not hasattr(minEnergy, '__len__'): - raise ValueError('minEnergy must be array-like if maxEnergy is') - if len(maxEnergy) != len(minEnergy): - raise ValueError('maxEnergy and minEnergy must be same length') - multipleEnergyRanges = True - lmin = [] - lmax = [] - for i in range(len(maxEnergy)): - lmin.append(ct.keVToAngstrom(maxEnergy[i])) - lmax.append(ct.keVToAngstrom(minEnergy[i])) - else: - lmin = ct.keVToAngstrom(maxEnergy) - lmax = ct.keVToAngstrom(minEnergy) - - # parse grain parameters kwarg - if grain_params is None: - grain_params = np.atleast_2d(np.hstack([np.zeros(6), ct.identity_6x1])) - n_grains = len(grain_params) - - # sample rotation - if rmat_s is None: - rmat_s = ct.identity_3x3 - - # dummy translation vector... make input - if tvec_s is None: - tvec_s = ct.zeros_3 - - # beam vector - if beam_vec is None: - beam_vec = ct.beam_vec - - # ========================================================================= - # LOOP OVER GRAINS - # ========================================================================= - - # pre-allocate output arrays - xy_det = np.nan * np.ones((n_grains, nhkls_tot, 2)) - hkls_in = np.nan * np.ones((n_grains, 3, nhkls_tot)) - angles = np.nan * np.ones((n_grains, nhkls_tot, 2)) - dspacing = np.nan * np.ones((n_grains, nhkls_tot)) - energy = np.nan * np.ones((n_grains, nhkls_tot)) - for iG, gp in enumerate(grain_params): - rmat_c = make_rmat_of_expmap(gp[:3]) - tvec_c = gp[3:6].reshape(3, 1) - vInv_s = mutil.vecMVToSymm(gp[6:].reshape(6, 1)) - - # stretch them: V^(-1) * R * Gc - gvec_s_str = np.dot(vInv_s, np.dot(rmat_c, gvec_c)) - ghat_c_str = mutil.unitVector(np.dot(rmat_c.T, gvec_s_str)) - - # project - dpts = gvec_to_xy( - ghat_c_str.T, - self.rmat, - rmat_s, - rmat_c, - self.tvec, - tvec_s, - tvec_c, - beam_vec=beam_vec, - ) - - # check intersections with detector plane - canIntersect = ~np.isnan(dpts[:, 0]) - npts_in = sum(canIntersect) - - if np.any(canIntersect): - dpts = dpts[canIntersect, :].reshape(npts_in, 2) - dhkl = hkls[:, canIntersect].reshape(3, npts_in) - - rmat_b = make_beam_rmat(beam_vec, ct.eta_vec) - # back to angles - tth_eta, gvec_l = xy_to_gvec( - dpts, - self.rmat, - rmat_s, - self.tvec, - tvec_s, - tvec_c, - rmat_b=rmat_b, - ) - tth_eta = np.vstack(tth_eta).T - - # warp measured points - if self.distortion is not None: - dpts = self.distortion.apply_inverse(dpts) - - # plane spacings and energies - dsp = 1.0 / mutil.rowNorm(gvec_s_str[:, canIntersect].T) - wlen = 2 * dsp * np.sin(0.5 * tth_eta[:, 0]) - - # clip to detector panel - _, on_panel = self.clip_to_panel(dpts, buffer_edges=True) - - if multipleEnergyRanges: - validEnergy = np.zeros(len(wlen), dtype=bool) - for i in range(len(lmin)): - in_energy_range = np.logical_and( - wlen >= lmin[i], wlen <= lmax[i] - ) - validEnergy = validEnergy | in_energy_range - else: - validEnergy = np.logical_and(wlen >= lmin, wlen <= lmax) - - # index for valid reflections - keepers = np.where(np.logical_and(on_panel, validEnergy))[0] - - # assign output arrays - xy_det[iG][keepers, :] = dpts[keepers, :] - hkls_in[iG][:, keepers] = dhkl[:, keepers] - angles[iG][keepers, :] = tth_eta[keepers, :] - dspacing[iG, keepers] = dsp[keepers] - energy[iG, keepers] = ct.keVToAngstrom(wlen[keepers]) - return xy_det, hkls_in, angles, dspacing, energy - @staticmethod def update_memoization_sizes(all_panels): funcs = [ diff --git a/hexrd/core/instrument/hedm_instrument.py b/hexrd/core/instrument/hedm_instrument.py index 8b5c54d9f..1a8c3fdc3 100644 --- a/hexrd/core/instrument/hedm_instrument.py +++ b/hexrd/core/instrument/hedm_instrument.py @@ -1547,52 +1547,6 @@ def simulate_powder_pattern( return img_dict - def simulate_laue_pattern( - self, - crystal_data, - minEnergy=5.0, - maxEnergy=35.0, - rmat_s=None, - grain_params=None, - ): - """ - Simulate Laue diffraction over the instrument. - - Parameters - ---------- - crystal_data : TYPE - DESCRIPTION. - minEnergy : TYPE, optional - DESCRIPTION. The default is 5.. - maxEnergy : TYPE, optional - DESCRIPTION. The default is 35.. - rmat_s : TYPE, optional - DESCRIPTION. The default is None. - grain_params : TYPE, optional - DESCRIPTION. The default is None. - - Returns - ------- - results : TYPE - DESCRIPTION. - - xy_det, hkls_in, angles, dspacing, energy - - TODO: revisit output; dict, or concatenated list? - """ - results = dict.fromkeys(self.detectors) - for det_key, panel in self.detectors.items(): - results[det_key] = panel.simulate_laue_pattern( - crystal_data, - minEnergy=minEnergy, - maxEnergy=maxEnergy, - rmat_s=rmat_s, - tvec_s=self.tvec, - grain_params=grain_params, - beam_vec=self.beam_vector, - ) - return results - def simulate_rotation_series( self, plane_data, diff --git a/hexrd/hedm/xrdutil/__init__.py b/hexrd/hedm/xrdutil/__init__.py index 3ad9db8e3..28c6d60b4 100644 --- a/hexrd/hedm/xrdutil/__init__.py +++ b/hexrd/hedm/xrdutil/__init__.py @@ -4,4 +4,3 @@ # TODO: Fully separate out the utils.py scripts from hexrd.hed.xrdutil.utils import * -from hexrd.laue.xrdutil.utils import * diff --git a/hexrd/laue/fitting/calibration/laue.py b/hexrd/laue/fitting/calibration/laue.py index 20e61c2f7..7b4ad47a7 100644 --- a/hexrd/laue/fitting/calibration/laue.py +++ b/hexrd/laue/fitting/calibration/laue.py @@ -15,6 +15,7 @@ from hexrd.core.rotations import angleAxisOfRotMat, RotMatEuler from hexrd.core.transforms import xfcapi from hexrd.core.utils.hkl import hkl_to_str, str_to_hkl +from hexrd.laue.simulation import simulate_laue_pattern_on_panel, simulate_laue_pattern_on_instrument # TODO: Resolve extra-workflow-dependency from ....core.fitting.calibration.calibrator import Calibrator @@ -183,7 +184,8 @@ def _autopick_points( # run simulation # ???: could we get this from overlays? - laue_sim = self.instr.simulate_laue_pattern( + laue_sim = simulate_laue_pattern_on_instrument( + self.instr, self.plane_data, minEnergy=self.energy_cutoffs[0], maxEnergy=self.energy_cutoffs[1], @@ -516,7 +518,8 @@ def sxcal_obj_func( for det_key, panel in instr.detectors.items(): # Simulate Laue pattern: # returns xy_det, hkls_in, angles, dspacing, energy - sim_results = panel.simulate_laue_pattern( + sim_results = simulate_laue_pattern_on_panel( + panel, [hkls_idx[det_key], bmat], minEnergy=energy_cutoffs[0], maxEnergy=energy_cutoffs[1], diff --git a/hexrd/laue/simulation.py b/hexrd/laue/simulation.py new file mode 100644 index 000000000..c80c36846 --- /dev/null +++ b/hexrd/laue/simulation.py @@ -0,0 +1,218 @@ +from __future__ import annotations + +import numpy as np +from numpy.typing import NDArray +from typing import TYPE_CHECKING, Optional, Sequence + +from hexrd.core.instrument.detector import Detector +from hexrd.core.material.crystallography import PlaneData + +if TYPE_CHECKING: + from hexrd.core.instrument.hedm_instrument import HEDMInstrument +import hexrd.core.constants as ct +from hexrd.core.transforms.xfcapi import gvec_to_xy, make_rmat_of_expmap, make_beam_rmat, xy_to_gvec +from hexrd.core import matrixutil as mutil + +def simulate_laue_pattern_on_panel( + detector: Detector, + crystal_data: PlaneData | Sequence[NDArray[np.float64]], + minEnergy: float | list[float] = 5.0, + maxEnergy: float | list[float] = 35.0, + rmat_s=None, + tvec_s=None, + grain_params=None, + beam_vec=None, +): + """ """ + if isinstance(crystal_data, PlaneData): + plane_data = crystal_data + + # grab the expanded list of hkls from plane_data + hkls = np.hstack(plane_data.getSymHKLs()) + + # and the unit plane normals (G-vectors) in CRYSTAL FRAME + gvec_c = np.dot(plane_data.latVecOps['B'], hkls) + + # Filter out g-vectors going in the wrong direction. `gvec_to_xy()` used + # to do this, but not anymore. + to_keep = np.dot(gvec_c.T, detector.bvec) <= 0 + + hkls = hkls[:, to_keep] + gvec_c = gvec_c[:, to_keep] + elif len(crystal_data) == 2: + # !!! should clean this up + hkls = np.array(crystal_data[0]) + bmat = crystal_data[1] + gvec_c = np.dot(bmat, hkls) + else: + raise RuntimeError(f'argument list not understood: {crystal_data=}') + nhkls_tot: int = hkls.shape[1] + + # parse energy ranges + # TODO: allow for spectrum parsing + multipleEnergyRanges = False + lmin: list[float] = [] + lmax: list[float] = [] + + if isinstance(maxEnergy, (list, tuple, np.ndarray)): + if not isinstance(minEnergy, (list, tuple, np.ndarray)): + raise TypeError('minEnergy must be array-like if maxEnergy is') + if len(maxEnergy) != len(minEnergy): + raise ValueError('maxEnergy and minEnergy must be same length') + multipleEnergyRanges = True + + for max_energy, min_energy in zip(maxEnergy, minEnergy): + lmin.append(ct.keVToAngstrom(max_energy)) + lmax.append(ct.keVToAngstrom(min_energy)) + else: + lmin = ct.keVToAngstrom(maxEnergy) + lmax = ct.keVToAngstrom(minEnergy) + + # parse grain parameters kwarg + if grain_params is None: + grain_params = np.atleast_2d(np.hstack([np.zeros(6), ct.identity_6x1])) + n_grains = len(grain_params) + + # sample rotation + if rmat_s is None: + rmat_s = ct.identity_3x3 + + # dummy translation vector... make input + if tvec_s is None: + tvec_s = ct.zeros_3 + + # beam vector + if beam_vec is None: + beam_vec = ct.beam_vec + + # ========================================================================= + # LOOP OVER GRAINS + # ========================================================================= + + # pre-allocate output arrays + xy_det = np.full((n_grains, nhkls_tot, 2), np.nan) + hkls_in = np.full((n_grains, 3, nhkls_tot), np.nan) + angles = np.full((n_grains, nhkls_tot, 2), np.nan) + dspacing = np.full((n_grains, nhkls_tot), np.nan) + energy = np.full((n_grains, nhkls_tot), np.nan) + + for iG, gp in enumerate(grain_params): + rmat_c = make_rmat_of_expmap(gp[:3]) + tvec_c = gp[3:6].reshape(3, 1) + vInv_s = mutil.vecMVToSymm(gp[6:].reshape(6, 1)) + + # stretch them: V^(-1) * R * Gc + gvec_s_str = np.dot(vInv_s, np.dot(rmat_c, gvec_c)) + ghat_c_str = mutil.unitVector(np.dot(rmat_c.T, gvec_s_str)) + + # project + dpts = gvec_to_xy( + ghat_c_str.T, + detector.rmat, + rmat_s, + rmat_c, + detector.tvec, + tvec_s, + tvec_c, + beam_vec=beam_vec, + ) + + # check intersections with detector plane + canIntersect = ~np.isnan(dpts[:, 0]) + npts_in = sum(canIntersect) + + if np.any(canIntersect): + dpts = dpts[canIntersect, :].reshape(npts_in, 2) + dhkl = hkls[:, canIntersect].reshape(3, npts_in) + + rmat_b = make_beam_rmat(beam_vec, ct.eta_vec) + # back to angles + tth_eta, gvec_l = xy_to_gvec( + dpts, + detector.rmat, + rmat_s, + detector.tvec, + tvec_s, + tvec_c, + rmat_b=rmat_b, + ) + tth_eta = np.vstack(tth_eta).T + + # warp measured points + if detector.distortion is not None: + dpts = detector.distortion.apply_inverse(dpts) + + # plane spacings and energies + dsp = 1.0 / mutil.rowNorm(gvec_s_str[:, canIntersect].T) + wlen = 2 * dsp * np.sin(0.5 * tth_eta[:, 0]) + + # clip to detector panel + _, on_panel = detector.clip_to_panel(dpts, buffer_edges=True) + + if multipleEnergyRanges: + validEnergy = np.zeros(len(wlen), dtype=bool) + for l_min, l_max in zip(lmin, lmax): + in_energy_range = np.logical_and( + wlen >= l_min, wlen <= l_max + ) + validEnergy = validEnergy | in_energy_range + else: + validEnergy = np.logical_and(wlen >= lmin, wlen <= lmax) + + # index for valid reflections + keepers = np.where(np.logical_and(on_panel, validEnergy))[0] + + # assign output arrays + xy_det[iG][keepers, :] = dpts[keepers, :] + hkls_in[iG][:, keepers] = dhkl[:, keepers] + angles[iG][keepers, :] = tth_eta[keepers, :] + dspacing[iG, keepers] = dsp[keepers] + energy[iG, keepers] = ct.keVToAngstrom(wlen[keepers]) + return xy_det, hkls_in, angles, dspacing, energy + +def simulate_laue_pattern_on_instrument( + instrument: HEDMInstrument, + crystal_data: PlaneData | Sequence[NDArray[np.float64]], + minEnergy: float = 5.0, + maxEnergy: float = 35.0, + rmat_s: Optional[NDArray[np.float64]] = None, + grain_params: Optional[NDArray[np.float64]] = None, +): + """ + Simulate Laue diffraction over the instrument. + + Parameters + ---------- + crystal_data : TYPE + DESCRIPTION. + minEnergy : TYPE, optional + DESCRIPTION. The default is 5.. + maxEnergy : TYPE, optional + DESCRIPTION. The default is 35.. + rmat_s : TYPE, optional + DESCRIPTION. The default is None. + grain_params : TYPE, optional + DESCRIPTION. The default is None. + + Returns + ------- + results : TYPE + DESCRIPTION. + + xy_det, hkls_in, angles, dspacing, energy + + TODO: revisit output; dict, or concatenated list? + """ + results = dict.fromkeys(instrument.detectors) + for det_key, panel in instrument.detectors.items(): + results[det_key] = simulate_laue_pattern_on_panel( + panel, + crystal_data, + minEnergy=minEnergy, + maxEnergy=maxEnergy, + rmat_s=rmat_s, + tvec_s=instrument.tvec, + grain_params=grain_params, + beam_vec=instrument.beam_vector, + ) + return results diff --git a/hexrd/laue/xrdutil/utils.py b/hexrd/laue/xrdutil/utils.py deleted file mode 100644 index 5c13a32e6..000000000 --- a/hexrd/laue/xrdutil/utils.py +++ /dev/null @@ -1,199 +0,0 @@ -#! /usr/bin/env python3 -# ============================================================ -# Copyright (c) 2012, Lawrence Livermore National Security, LLC. -# Produced at the Lawrence Livermore National Laboratory. -# Written by Joel Bernier and others. -# LLNL-CODE-529294. -# All rights reserved. -# -# This file is part of HEXRD. For details on dowloading the source, -# see the file COPYING. -# -# Please also see the file LICENSE. -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License (as published by the Free -# Software Foundation) version 2.1 dated February 1999. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY -# or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this program (see file LICENSE); if not, write to -# the Free Software Foundation, Inc., 59 Temple Place, Suite 330, -# Boston, MA 02111-1307 USA or visit . -# ============================================================ - - -import numpy as np - -from hexrd.core import constants -from hexrd.core import matrixutil as mutil - -from hexrd.core.material.crystallography import processWavelength - -from hexrd.core.transforms import xfcapi - -from hexrd.core.deprecation import deprecated - - -simlp = 'hexrd.core.instrument.hedm_instrument.HEDMInstrument.simulate_laue_pattern' - -# ============================================================================= -# PARAMETERS -# ============================================================================= - -distortion_key = 'distortion' - -d2r = piby180 = constants.d2r -r2d = constants.r2d - -epsf = constants.epsf # ~2.2e-16 -ten_epsf = 10 * epsf # ~2.2e-15 -sqrt_epsf = constants.sqrt_epsf # ~1.5e-8 - -bHat_l_DFLT = constants.beam_vec.flatten() -eHat_l_DFLT = constants.eta_vec.flatten() - -nans_1x2 = np.nan * np.ones((1, 2)) - -validateAngleRanges = xfcapi.validate_angle_ranges - - -@deprecated(new_func=simlp, removal_date='2026-01-01') -def simulateLauePattern( - hkls, - bMat, - rmat_d, - tvec_d, - panel_dims, - panel_buffer=5, - minEnergy=8, - maxEnergy=24, - rmat_s=np.eye(3), - grain_params=None, - distortion=None, - beamVec=None, -): - - if beamVec is None: - beamVec = constants.beam_vec - - # parse energy ranges - multipleEnergyRanges = False - if hasattr(maxEnergy, '__len__'): - assert len(maxEnergy) == len( - minEnergy - ), 'energy cutoff ranges must have the same length' - multipleEnergyRanges = True - lmin = [processWavelength(e) for e in maxEnergy] - lmax = [processWavelength(e) for e in minEnergy] - else: - lmin = processWavelength(maxEnergy) - lmax = processWavelength(minEnergy) - - # process crystal rmats and inverse stretches - if grain_params is None: - grain_params = np.atleast_2d( - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0] - ) - - n_grains = len(grain_params) - - # dummy translation vector... make input - tvec_s = np.zeros((3, 1)) - - # number of hkls - nhkls_tot = hkls.shape[1] - - # unit G-vectors in crystal frame - ghat_c = mutil.unitVector(np.dot(bMat, hkls)) - - # pre-allocate output arrays - xy_det = np.nan * np.ones((n_grains, nhkls_tot, 2)) - hkls_in = np.nan * np.ones((n_grains, 3, nhkls_tot)) - angles = np.nan * np.ones((n_grains, nhkls_tot, 2)) - dspacing = np.nan * np.ones((n_grains, nhkls_tot)) - energy = np.nan * np.ones((n_grains, nhkls_tot)) - - """ - LOOP OVER GRAINS - """ - - for iG, gp in enumerate(grain_params): - rmat_c = xfcapi.make_rmat_of_expmap(gp[:3]) - tvec_c = gp[3:6].reshape(3, 1) - vInv_s = mutil.vecMVToSymm(gp[6:].reshape(6, 1)) - - # stretch them: V^(-1) * R * Gc - ghat_s_str = mutil.unitVector(np.dot(vInv_s, np.dot(rmat_c, ghat_c))) - ghat_c_str = np.dot(rmat_c.T, ghat_s_str) - - # project - dpts = xfcapi.gvec_to_xy( - ghat_c_str.T, - rmat_d, - rmat_s, - rmat_c, - tvec_d, - tvec_s, - tvec_c, - beam_vec=beamVec, - ).T - - # check intersections with detector plane - canIntersect = ~np.isnan(dpts[0, :]) - npts_in = sum(canIntersect) - - if np.any(canIntersect): - dpts = dpts[:, canIntersect].reshape(2, npts_in) - dhkl = hkls[:, canIntersect].reshape(3, npts_in) - - rmat_b = xfcapi.make_beam_rmat(beamVec, constants.eta_vec) - - # back to angles - tth_eta, gvec_l = xfcapi.xy_to_gvec( - dpts.T, rmat_d, rmat_s, tvec_d, tvec_s, tvec_c, rmat_b=rmat_b - ) - tth_eta = np.vstack(tth_eta).T - - # warp measured points - if distortion is not None: - dpts = distortion.apply_inverse(dpts) - - # plane spacings and energies - dsp = 1.0 / mutil.columnNorm(np.dot(bMat, dhkl)) - wlen = 2 * dsp * np.sin(0.5 * tth_eta[:, 0]) - - # find on spatial extent of detector - xTest = np.logical_and( - dpts[0, :] >= -0.5 * panel_dims[1] + panel_buffer, - dpts[0, :] <= 0.5 * panel_dims[1] - panel_buffer, - ) - yTest = np.logical_and( - dpts[1, :] >= -0.5 * panel_dims[0] + panel_buffer, - dpts[1, :] <= 0.5 * panel_dims[0] - panel_buffer, - ) - - onDetector = np.logical_and(xTest, yTest) - if multipleEnergyRanges: - validEnergy = np.zeros(len(wlen), dtype=bool) - for i in range(len(lmin)): - validEnergy = validEnergy | np.logical_and( - wlen >= lmin[i], wlen <= lmax[i] - ) - else: - validEnergy = np.logical_and(wlen >= lmin, wlen <= lmax) - - # index for valid reflections - keepers = np.where(np.logical_and(onDetector, validEnergy))[0] - - # assign output arrays - xy_det[iG][keepers, :] = dpts[:, keepers].T - hkls_in[iG][:, keepers] = dhkl[:, keepers] - angles[iG][keepers, :] = tth_eta[keepers, :] - dspacing[iG, keepers] = dsp[keepers] - energy[iG, keepers] = processWavelength(wlen[keepers]) - return xy_det, hkls_in, angles, dspacing, energy diff --git a/tests/core/fitting/calibration/test_fitting_calibration_laue.py b/tests/core/fitting/calibration/test_fitting_calibration_laue.py index 46ddda7bb..3161e690e 100644 --- a/tests/core/fitting/calibration/test_fitting_calibration_laue.py +++ b/tests/core/fitting/calibration/test_fitting_calibration_laue.py @@ -97,7 +97,9 @@ def test_laue_attributes_and_setters(laue_calibrator): @patch('hexrd.core.fitting.calibration.laue.xrdutil') @patch('hexrd.core.fitting.calibration.laue.xfcapi') @patch('hexrd.core.fitting.calibration.laue.switch_xray_source') +@patch('hexrd.core.fitting.calibration.laue.simulate_laue_pattern_on_instrument') def test_autopick_points( + mock_sim_instr, mock_switch, mock_xfc, mock_xrd, @@ -123,7 +125,7 @@ def test_autopick_points( "lsq", ([100, 1, 0, 1, 5.0, 5.0, 0, 0, 0], 1) ) - laue_calibrator.instr.simulate_laue_pattern.return_value = { + mock_sim_instr.return_value = { 'det1': ( np.array([[[100, 100]]]), np.array([[[1], [1], [1]]]), @@ -176,10 +178,11 @@ def test_autopick_points( # --- Objective Function & Methods --- +@patch('hexrd.core.fitting.calibration.laue.simulate_laue_pattern_on_panel') @patch('hexrd.core.fitting.calibration.laue.switch_xray_source') @patch('hexrd.core.fitting.calibration.laue.sxcal_obj_func') def test_laue_methods_and_obj_func( - mock_obj_func, mock_switch, laue_calibrator + mock_obj_func, mock_switch, mock_sim, laue_calibrator ): laue_calibrator.data_dict = { 'pick_xys': {'det1': [[10, 10], [np.nan, np.nan]]}, @@ -196,12 +199,12 @@ def test_laue_methods_and_obj_func( assert mock_obj_func.call_args[1]['sim_only'] is True mock_instr = laue_calibrator.instr - mock_instr.detectors['det1'].simulate_laue_pattern.return_value = ( - [np.array([[100, 100]])], - [np.zeros((1, 3))], - [np.array([[10, 20]])], - None, - None, + mock_sim.return_value = ( + np.array([[[100, 100]]]), # xy_det: (n_grains, nhkls, 2) + np.zeros((1, 3, 1)), # hkls_in: (n_grains, 3, nhkls) + np.array([[[10, 20]]]), # angles: (n_grains, nhkls, 2) + np.array([[1.0]]), # dspacing: (n_grains, nhkls) + np.array([[20.0]]), # energy: (n_grains, nhkls) ) res_sim = sxcal_obj_func( diff --git a/tests/core/instrument/test_detector.py b/tests/core/instrument/test_detector.py index 2ae142507..6405af5bf 100644 --- a/tests/core/instrument/test_detector.py +++ b/tests/core/instrument/test_detector.py @@ -10,6 +10,7 @@ _interpolate_bilinear, _interpolate_bilinear_in_place, ) +from hexrd.laue.simulation import simulate_laue_pattern_on_panel from hexrd.core.instrument.physics_package import AbstractPhysicsPackage from hexrd.core.material import crystallography @@ -570,34 +571,39 @@ def test_simulate_laue_pattern(base_detector, mock_distortion_registry): pd.getSymHKLs.return_value = [np.array([[1], [0], [0]])] pd.latVecOps = {'B': np.eye(3)} - with patch('hexrd.core.instrument.detector.xy_to_gvec') as m_xy2g: + with patch('hexrd.laue.simulation.xy_to_gvec') as m_xy2g: m_xy2g.return_value = ([np.array([0.1, 0.0])], np.array([[0, 0, 1]])) with patch( - 'hexrd.core.instrument.detector.gvec_to_xy', + 'hexrd.laue.simulation.gvec_to_xy', return_value=np.array([[0.0, 0.0]]), ): with np.errstate(divide='ignore', invalid='ignore'): - base_detector.simulate_laue_pattern( - pd, grain_params=[np.zeros(12)] + simulate_laue_pattern_on_panel( + base_detector, pd, grain_params=[np.zeros(12)] ) with pytest.raises( - ValueError, + TypeError, match="minEnergy must be array-like if maxEnergy is", ): - base_detector.simulate_laue_pattern(pd, maxEnergy=[1, 2]) + simulate_laue_pattern_on_panel( + base_detector, pd, maxEnergy=[1, 2] + ) with pytest.raises( ValueError, match="maxEnergy and minEnergy must be same length", ): - base_detector.simulate_laue_pattern( - pd, minEnergy=[1], maxEnergy=[1, 2] + simulate_laue_pattern_on_panel( + base_detector, pd, minEnergy=[1], maxEnergy=[1, 2] ) - base_detector.simulate_laue_pattern( - pd, minEnergy=[5.0, 10.0], maxEnergy=[15.0, 20.0] + simulate_laue_pattern_on_panel( + base_detector, + pd, + minEnergy=[5.0, 10.0], + maxEnergy=[15.0, 20.0], ) @@ -606,21 +612,21 @@ def test_simulate_laue_pattern_legacy_input(base_detector): bmat = np.eye(3) with patch( - 'hexrd.core.instrument.detector.gvec_to_xy', + 'hexrd.laue.simulation.gvec_to_xy', return_value=np.array([[0.0, 0.0]]), ): - with patch('hexrd.core.instrument.detector.xy_to_gvec') as m_xy2g: + with patch('hexrd.laue.simulation.xy_to_gvec') as m_xy2g: m_xy2g.return_value = ( [np.array([0.1, 0.0])], np.array([[0, 0, 1]]), ) with np.errstate(divide='ignore', invalid='ignore'): - base_detector.simulate_laue_pattern((hkls, bmat)) + simulate_laue_pattern_on_panel(base_detector, (hkls, bmat)) def test_simulate_laue_pattern_error(base_detector): with pytest.raises(RuntimeError): - base_detector.simulate_laue_pattern("invalid_input") + simulate_laue_pattern_on_panel(base_detector, "invalid_input") def test_polarization_factor(base_detector): diff --git a/tests/test_laue.py b/tests/test_laue.py index 26f26d6f2..1c2e2c18c 100644 --- a/tests/test_laue.py +++ b/tests/test_laue.py @@ -7,6 +7,7 @@ from hexrd.core.material.material import load_materials_hdf5, Material from hexrd.core.instrument.hedm_instrument import HEDMInstrument +from hexrd.laue.simulation import simulate_laue_pattern_on_instrument @pytest.fixture @@ -54,8 +55,8 @@ def test_simulate_laue_spots( # Disable all exclusions on the plane data plane_data.exclusions = None - sim_data = instr.simulate_laue_pattern( - plane_data, minEnergy=5, maxEnergy=35, grain_params=[lif_grain_params] + sim_data = simulate_laue_pattern_on_instrument( + instr, plane_data, minEnergy=5, maxEnergy=35, grain_params=[lif_grain_params] ) # A few manual expected results