From c8e53a19350e4571f3249bb863838c64d3547fb9 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 4 Mar 2024 12:24:20 -0600 Subject: [PATCH] Add hedm_intensity to PlaneData This is computed nearly the same way as the powder intensity, except that it does not include multiplicity. This PR also moves the lorentz and polarization functions to the crystallography module, where they are now used for computing the powder and HEDM intensities. Signed-off-by: Patrick Avery --- hexrd/instrument/detector.py | 71 +++++++------------------------ hexrd/material/crystallography.py | 68 +++++++++++++++++++++++++++-- 2 files changed, 81 insertions(+), 58 deletions(-) diff --git a/hexrd/instrument/detector.py b/hexrd/instrument/detector.py index 6cfd8d834..a05fc6c9c 100644 --- a/hexrd/instrument/detector.py +++ b/hexrd/instrument/detector.py @@ -9,6 +9,7 @@ from hexrd import matrixutil as mutil from hexrd import xrdutil +from hexrd.material import crystallography from hexrd.material.crystallography import PlaneData from hexrd.transforms.xfcapi import ( @@ -38,6 +39,12 @@ beam_energy_DFLT = 65.351 +# Memoize these, so each detector can avoid re-computing if nothing +# has changed. +_lorentz_factor = memoize(crystallography.lorentz_factor) +_polarization_factor = memoize(crystallography.polarization_factor) + + class Detector: """ Base class for 2D detectors with functions and properties @@ -551,9 +558,15 @@ def polarization_factor(self, f_hor, f_vert, unpolarized=False): raise RuntimeError(msg) tth, eta = self.pixel_angles() - args = (tth, eta, f_hor, f_vert, unpolarized) + kwargs = { + 'tth': tth, + 'eta': eta, + 'f_hor': f_hor, + 'f_vert': f_vert, + 'unpolarized': unpolarized, + } - return _polarization_factor(*args) + return _polarization_factor(**kwargs) def lorentz_factor(self): """ @@ -575,9 +588,7 @@ def lorentz_factor(self): corresponding pixel """ tth, eta = self.pixel_angles() - args = (tth,) - - return _lorentz_factor(*args) + return _lorentz_factor(tth) def config_dict(self, chi=0, tvec=ct.zeros_3, beam_energy=beam_energy_DFLT, beam_vector=ct.beam_vec, @@ -1548,56 +1559,6 @@ def _col_edge_vec(cols, pixel_size_col): return pixel_size_col*(np.arange(cols+1)-0.5*cols) -@memoize -def _polarization_factor(tth, eta, f_hor, f_vert, unpolarized): - """ - 06/14/2021 SS adding lorentz polarization factor computation - to the detector so that it can be compenstated for in the - intensity correction - - 05/26/2022 decoupling lorentz factor from polarization factor - - parameters: tth two theta of every pixel in radians - eta azimuthal angle of every pixel - f_hor fraction of horizontal polarization - (~1 for XFELs) - f_vert fraction of vertical polarization - (~0 for XFELs) - notice f_hor + f_vert = 1 - """ - - ctth2 = np.cos(tth)**2 - seta2 = np.sin(eta)**2 - ceta2 = np.cos(eta)**2 - - if unpolarized: - P = (1. + ctth2)/2. - else: - P = f_hor*(seta2 + ceta2*ctth2) + f_vert*(ceta2 + seta2*ctth2) - - return P - - -@memoize -def _lorentz_factor(tth): - """ - 05/26/2022 SS adding lorentz factor computation - to the detector so that it can be compenstated for in the - intensity correction - - parameters: tth two theta of every pixel in radians - """ - - theta = 0.5*tth - - cth = np.cos(theta) - sth2 = np.sin(theta)**2 - - L = 1./(4.0*cth*sth2) - - return L - - # FIXME find a better place for this, and maybe include loop over pixels if ct.USE_NUMBA: @numba.njit(nogil=True, cache=True) diff --git a/hexrd/material/crystallography.py b/hexrd/material/crystallography.py index d3d5c63c4..2ab5b1738 100644 --- a/hexrd/material/crystallography.py +++ b/hexrd/material/crystallography.py @@ -943,12 +943,14 @@ def invalidate_structure_factor(self, unitcell): # option to just invalidate it, while providing a unit cell, so that # it can be lazily computed from the unit cell. self.__structFact = None + self._hedm_intensity = None self._powder_intensity = None self.__unitcell = unitcell def _compute_sf_if_needed(self): any_invalid = ( self.__structFact is None or + self._hedm_intensity is None or self._powder_intensity is None ) if any_invalid and self.__unitcell is not None: @@ -965,11 +967,18 @@ def set_structFact(self, structFact): self.__structFact = structFact multiplicity = self.getMultiplicity(allHKLs=True) tth = self.getTTh(allHKLs=True) - lp = (1 + np.cos(tth)**2)/np.cos(0.5*tth)/np.sin(0.5*tth)**2/2.0 - powderI = structFact*multiplicity*lp - powderI = 100.0*powderI/np.nanmax(powderI) + hedm_intensity = ( + structFact * lorentz_factor(tth) * polarization_factor(tth) + ) + + powderI = hedm_intensity * multiplicity + + # Now scale them + hedm_intensity = 100.0 * hedm_intensity / np.nanmax(hedm_intensity) + powderI = 100.0 * powderI / np.nanmax(powderI) + self._hedm_intensity = hedm_intensity self._powder_intensity = powderI structFact = property(get_structFact, set_structFact, None) @@ -979,6 +988,11 @@ def powder_intensity(self): self._compute_sf_if_needed() return self._powder_intensity[~self.exclusions] + @property + def hedm_intensity(self): + self._compute_sf_if_needed() + return self._hedm_intensity[~self.exclusions] + @staticmethod def makePlaneData(hkls, lparms, qsym, symmGroup, strainMag, wavelength): """ @@ -1947,3 +1961,51 @@ def RetrieveAtomicFormFactor(elecNum, magG, sinThOverLamdaList, ffDataList): ) return ff + + +def lorentz_factor(tth): + """ + 05/26/2022 SS adding lorentz factor computation + to the detector so that it can be compenstated for in the + intensity correction + + parameters: tth two theta of every pixel in radians + """ + + theta = 0.5*tth + + cth = np.cos(theta) + sth2 = np.sin(theta)**2 + + L = 1./(4.0*cth*sth2) + + return L + + +def polarization_factor(tth, unpolarized=True, eta=None, + f_hor=None, f_vert=None): + """ + 06/14/2021 SS adding lorentz polarization factor computation + to the detector so that it can be compenstated for in the + intensity correction + + 05/26/2022 decoupling lorentz factor from polarization factor + + parameters: tth two theta of every pixel in radians + if unpolarized is True, all subsequent arguments are optional + eta azimuthal angle of every pixel + f_hor fraction of horizontal polarization + (~1 for XFELs) + f_vert fraction of vertical polarization + (~0 for XFELs) + notice f_hor + f_vert = 1 + """ + + ctth2 = np.cos(tth)**2 + + if unpolarized: + return (1 + ctth2) / 2 + + seta2 = np.sin(eta)**2 + ceta2 = np.cos(eta)**2 + return f_hor*(seta2 + ceta2*ctth2) + f_vert*(ceta2 + seta2*ctth2)