Merge pull request #619 from HEXRD/hedm-intensity

Add hedm_intensity to PlaneData

.github/workflows/container_build.sh

@@ -16,6 +16,9 @@ wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
 chmod a+x Miniconda3-latest-Linux-x86_64.sh
 ./Miniconda3-latest-Linux-x86_64.sh -b
 
+# The base needs to have the same python version (3.11, right now)
+${HOME}/miniconda3/bin/conda install python=3.11 -y
+
 # Set up the hexrd channel
 ${HOME}/miniconda3/bin/conda create --override-channels -c conda-forge -y -n hexrd python=3.11
 ${HOME}/miniconda3/bin/activate hexrd

.github/workflows/package.yml

@@ -36,7 +36,7 @@ jobs:
         fetch-depth: 0
 
     - name: Install conda
-      uses: conda-incubator/setup-miniconda@v2
+      uses: conda-incubator/setup-miniconda@v3
       with:
         auto-update-conda: true
         python-version: 3.11

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)