Skip to content

Simulate laue pattern #911

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ZackAttack614 wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Simulate laue pattern #911

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fe2aeed
Move simulate_laue_pattern into the laue workflow
ZackAttack614 Mar 26, 2026
51f83f5
Delete old simulateLauePattern
ZackAttack614 Mar 26, 2026
c3cf46b
Add import for simulate_laue_pattern_on_instrument
ZackAttack614 Mar 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions hexrd/core/fitting/calibration/laue.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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],
Expand Down
155 changes: 0 additions & 155 deletions hexrd/core/instrument/detector.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = [
Expand Down
46 changes: 0 additions & 46 deletions hexrd/core/instrument/hedm_instrument.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
1 change: 0 additions & 1 deletion hexrd/hedm/xrdutil/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,4 +4,3 @@

# TODO: Fully separate out the utils.py scripts
from hexrd.hed.xrdutil.utils import *
from hexrd.laue.xrdutil.utils import *
7 changes: 5 additions & 2 deletions hexrd/laue/fitting/calibration/laue.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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],
Expand Down
Loading
Loading