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 <[email protected]>

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)

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)