feat: expose ReflectionProfileDoubleExponentialPseudoVoigt TOF profile

news/tof-profile.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* Expose the time-of-flight reflection profile ReflectionProfileDoubleExponentialPseudoVoigt.
+
+**Changed:**
+
+* PowderPatternDiffraction.SetProfile gives a ReflectionProfileDoubleExponentialPseudoVoigt the unit cell of the component's crystal.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/extensions/powderpatterndiffraction_ext.cpp

@@ -41,7 +41,14 @@ namespace
 void SetProfileCopy(PowderPatternDiffraction &pdiff,
                     const ReflectionProfile &profile)
 {
-    pdiff.SetProfile(profile.CreateCopy());
+    ReflectionProfile *copy = profile.CreateCopy();
+    // The TOF profile needs a unit cell to compute d_hkl, as in
+    // PowderPatternDiffraction::SetCrystal and XMLInput.
+    ReflectionProfileDoubleExponentialPseudoVoigt *tofprofile =
+        dynamic_cast<ReflectionProfileDoubleExponentialPseudoVoigt *>(copy);
+    if (tofprofile != 0 && pdiff.HasCrystal())
+        tofprofile->SetUnitCell(pdiff.GetCrystal());
+    pdiff.SetProfile(copy);
 }
 
 // PowderPatternDiffraction has no public accessor for its texture

src/extensions/reflectionprofile_ext.cpp

@@ -159,6 +159,31 @@ void wrap_reflectionprofile()
               bp::arg("scherrerLGmix") = 1),
              "Set the TCH pseudo-Voigt U, V, W, X, Y, Z, P and LGmix parameters.");
 
+    class_<ReflectionProfileDoubleExponentialPseudoVoigt,
+           bases<ReflectionProfile> >(
+        "ReflectionProfileDoubleExponentialPseudoVoigt", init<>())
+        .def(init<const UnitCell &>(bp::arg("cell"))
+             [with_custodian_and_ward<1, 2>()])
+        .def(init<const ReflectionProfileDoubleExponentialPseudoVoigt &>(
+            bp::arg("old"))[with_custodian_and_ward<1, 2>()])
+        .def("CreateCopy",
+             &ReflectionProfileDoubleExponentialPseudoVoigt::CreateCopy,
+             return_value_policy<manage_new_object,
+                                 with_custodian_and_ward_postcall<0, 1> >(),
+             "Return an independent copy of this profile.")
+        .def("SetProfilePar",
+             &ReflectionProfileDoubleExponentialPseudoVoigt::SetProfilePar,
+             (bp::arg("instrumentAlpha0"), bp::arg("instrumentAlpha1"),
+              bp::arg("instrumentBeta0"), bp::arg("instrumentBeta1"),
+              bp::arg("gaussianSigma0"), bp::arg("gaussianSigma1"),
+              bp::arg("gaussianSigma2"), bp::arg("lorentzianGamma0"),
+              bp::arg("lorentzianGamma1"), bp::arg("lorentzianGamma2")),
+             "Set the double-exponential pseudo-Voigt TOF profile parameters.")
+        .def("SetUnitCell",
+             &ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell,
+             bp::arg("cell"), with_custodian_and_ward<1, 2>(),
+             "Set the unit cell used to compute d_hkl for reflections.");
+
     class_<ReflectionProfilePseudoVoigtAnisotropic,
            bases<ReflectionProfile> >(
         "ReflectionProfilePseudoVoigtAnisotropic", init<>())

src/pyobjcryst/reflectionprofile.py

@@ -36,6 +36,11 @@
     for pdiff in powder_pattern.get_crystalline_components():
         pdiff.SetProfile(profile)
 
+For time-of-flight neutron data, use the double-exponential pseudo-Voigt
+profile ``ReflectionProfileDoubleExponentialPseudoVoigt``. It needs the
+unit cell to compute ``d_hkl``, which ``PowderPatternDiffraction.SetProfile``
+sets from the component's crystal.
+
 Example
 -------
 
@@ -69,11 +74,13 @@
     "ReflectionProfilePseudoVoigt",
     "ReflectionProfilePseudoVoigtTCH",
     "ReflectionProfilePseudoVoigtAnisotropic",
+    "ReflectionProfileDoubleExponentialPseudoVoigt",
     "ReflectionProfileType",
 ]
 
 from pyobjcryst._pyobjcryst import (
     ReflectionProfile,
+    ReflectionProfileDoubleExponentialPseudoVoigt,
     ReflectionProfilePseudoVoigt,
     ReflectionProfilePseudoVoigtAnisotropic,
     ReflectionProfilePseudoVoigtTCH,

tests/test_powderpattern.py

@@ -483,7 +483,7 @@ def test_add_phase(loadcifdata, phases, expected_par_names):
         texture.GetPar(i).GetName() for i in range(texture.GetNbPar())
     ]
     assert actual_type == expected_type
-    assert actual_phases == pytest.approx(expected_phases)
+    assert actual_phases == [pytest.approx(e) for e in expected_phases]
     assert actual_par_names == expected_par_names
 
 

tests/test_reflectionprofile.py

@@ -1,14 +1,18 @@
 """Unit tests for pyobjcryst.reflectionprofile bindings."""
 
+import gc
 import unittest
+import weakref
 
 import numpy as np
 import pytest
 
 from pyobjcryst.powderpattern import PowderPattern
+from pyobjcryst.radiation import RadiationType, WavelengthType
 from pyobjcryst.refinableobj import RefinableObj
 from pyobjcryst.reflectionprofile import (
     ReflectionProfile,
+    ReflectionProfileDoubleExponentialPseudoVoigt,
     ReflectionProfilePseudoVoigt,
     ReflectionProfilePseudoVoigtAnisotropic,
     ReflectionProfilePseudoVoigtTCH,
@@ -293,5 +297,140 @@ def test_xml_output(self):
         self.assertTrue(xml_from_stream.startswith("<ReflectionProfile"))
 
 
+def build_tof_pattern(crystal, profile):
+    """Return a TOF neutron powder pattern of crystal using profile."""
+    pattern = PowderPattern()
+    pattern.SetRadiationType(RadiationType.RAD_NEUTRON)
+    pattern.GetRadiation().SetWavelengthType(WavelengthType.WAVELENGTH_TOF)
+    # TOF decreases with sin(theta)/lambda, so x is in descending order.
+    x = np.linspace(120000, 20000, 2001)
+    pattern.SetPowderPatternX(x)
+    pattern.SetPowderPatternObs(np.ones_like(x))
+    diffraction = pattern.AddPowderPatternDiffraction(crystal)
+    diffraction.SetProfile(profile)
+    return pattern
+
+
+def test_set_profile_par_tof():
+    # Test that SetProfilePar assigns the ten TOF profile parameters
+    # Expected: Each refinable parameter has the value passed for it.
+    profile = ReflectionProfileDoubleExponentialPseudoVoigt()
+    profile.SetProfilePar(
+        instrumentAlpha0=0.01,
+        instrumentAlpha1=0.1,
+        instrumentBeta0=0.02,
+        instrumentBeta1=0.004,
+        gaussianSigma0=1.0,
+        gaussianSigma1=7.0,
+        gaussianSigma2=0.5,
+        lorentzianGamma0=0.1,
+        lorentzianGamma1=0.2,
+        lorentzianGamma2=0.4,
+    )
+    actual_values = {
+        profile.GetPar(i).GetName(): profile.GetPar(i).GetValue()
+        for i in range(profile.GetNbPar())
+    }
+    expected_values = {
+        "Alpha0": 0.01,
+        "Alpha1": 0.1,
+        "Beta0": 0.02,
+        "Beta1": 0.004,
+        "GaussianSigma0": 1.0,
+        "GaussianSigma1": 7.0,
+        "GaussianSigma2": 0.5,
+        "LorentzianGamma0": 0.1,
+        "LorentzianGamma1": 0.2,
+        "LorentzianGamma2": 0.4,
+    }
+    assert actual_values == pytest.approx(expected_values)
+
+
+def test_set_profile_tof_without_unit_cell(loadcifdata):
+    # Test installing a TOF profile constructed without a unit cell
+    # Expected: SetProfile gives the profile the crystal's unit cell, so
+    # the strongest peak is the Si (220) reflection at
+    # TOF = DIFC * d + DIFA * d^2.
+    crystal = loadcifdata("Si.cif")
+    pattern = build_tof_pattern(
+        crystal, ReflectionProfileDoubleExponentialPseudoVoigt()
+    )
+    calc = pattern.GetPowderPatternCalc()
+    x = pattern.GetPowderPatternX()
+    actual_peak_tof = x[np.argmax(calc)]
+    d_220 = crystal.a / np.sqrt(8)
+    expected_peak_tof = (
+        pattern.GetPar("DIFC").GetValue() * d_220
+        + pattern.GetPar("DIFA").GetValue() * d_220**2
+    )
+    tof_step = abs(x[1] - x[0])
+    assert actual_peak_tof == pytest.approx(expected_peak_tof, abs=tof_step)
+
+
+@pytest.mark.parametrize(
+    "parameter_name, new_value",
+    [
+        # Test that changing a parameter of the installed TOF profile
+        # recalculates the pattern
+        # C1: The rising-edge exponential coefficient Alpha1 is changed.
+        # Expected: The pattern equals one calculated with Alpha1 set
+        # before the profile was installed.
+        ("Alpha1", 0.2),
+        # C2: The falling-edge exponential coefficient Beta0 is changed.
+        # Expected: The pattern equals one calculated with Beta0 set
+        # before the profile was installed.
+        ("Beta0", 0.05),
+        # C3: The Gaussian width coefficient GaussianSigma1 is changed.
+        # Expected: The pattern equals one calculated with GaussianSigma1
+        # set before the profile was installed.
+        ("GaussianSigma1", 50.0),
+        # C4: The Lorentzian width coefficient LorentzianGamma2 is changed.
+        # Expected: The pattern equals one calculated with
+        # LorentzianGamma2 set before the profile was installed.
+        ("LorentzianGamma2", 2.0),
+    ],
+)
+def test_get_powder_pattern_calc_after_tof_profile_change(
+    loadcifdata, parameter_name, new_value
+):
+    crystal = loadcifdata("Si.cif")
+    pattern = build_tof_pattern(
+        crystal, ReflectionProfileDoubleExponentialPseudoVoigt()
+    )
+    pattern.GetPowderPatternCalc()
+    profile = pattern.GetPowderPatternComponent(0).GetProfile()
+    profile.GetPar(parameter_name).SetValue(new_value)
+    actual_calc = pattern.GetPowderPatternCalc()
+    expected_profile = ReflectionProfileDoubleExponentialPseudoVoigt(crystal)
+    expected_profile.GetPar(parameter_name).SetValue(new_value)
+    expected_calc = build_tof_pattern(
+        crystal, expected_profile
+    ).GetPowderPatternCalc()
+    assert np.allclose(actual_calc, expected_calc)
+
+
+@pytest.mark.parametrize("copy_method", ["CreateCopy", "copy_constructor"])
+def test_tof_profile_copy_keeps_original_alive(loadcifdata, copy_method):
+    # Test that a copy of a TOF profile keeps its source and unit cell alive
+    # Expected: The original profile and crystal survive until the copy is
+    # deleted, since both profiles reference the same unit cell.
+    crystal = loadcifdata("Si.cif")
+    original = ReflectionProfileDoubleExponentialPseudoVoigt(crystal)
+    if copy_method == "CreateCopy":
+        copied = original.CreateCopy()
+    else:
+        copied = ReflectionProfileDoubleExponentialPseudoVoigt(original)
+    original_ref = weakref.ref(original)
+    crystal_ref = weakref.ref(crystal)
+    del original, crystal
+    gc.collect()
+    assert original_ref() is not None
+    assert crystal_ref() is not None
+    del copied
+    gc.collect()
+    assert original_ref() is None
+    assert crystal_ref() is None
+
+
 if __name__ == "__main__":
     unittest.main()