Refs #22420 add cross terms to predict

Framework/Crystal/inc/MantidCrystal/PredictSatellitePeaks.h

@@ -67,12 +67,19 @@ private:
   void exec_peaks();
   Kernel::V3D getOffsetVector(const std::string &label);
   void predictOffsets(DataObjects::PeaksWorkspace_sptr Peaks,
-                      boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks,
+                      boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, int iVector,
                       Kernel::V3D offsets, int &maxOrder, Kernel::V3D &hkl,
                       Geometry::HKLFilterWavelength &lambdaFilter,
                       bool &includePeaksInRange, bool &includeOrderZero,
                       int &seqNum,
                       std::vector<std::vector<int>> &AlreadyDonePeaks);
+void predictOffsetsWithCrossTerms(
+    DataObjects::PeaksWorkspace_sptr Peaks,
+    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, Kernel::V3D offsets1,
+    Kernel::V3D offsets2, Kernel::V3D offsets3,
+    int &maxOrder, Kernel::V3D &hkl, Geometry::HKLFilterWavelength &lambdaFilter,
+    bool &includePeaksInRange, bool &includeOrderZero, int &seqNum,
+    std::vector<std::vector<int>> &AlreadyDonePeaks);
 };
 
 } // namespace Crystal

Framework/Crystal/src/PredictSatellitePeaks.cpp

@@ -53,6 +53,10 @@ void PredictSatellitePeaks::init() {
       make_unique<PropertyWithValue<int>>("MaxOrder", 0, Direction::Input),
       "Maximum order to apply ModVectors. Default = 0");
 
+  declareProperty(make_unique<PropertyWithValue<bool>>("CrossTerms", false,
+                                                       Direction::Input),
+                  "Include cross terms (false)");
+
   declareProperty(
       "IncludeIntegerHKL", true,
       "If false order 0 peaks are not included in workspace (integer HKL)");
@@ -113,6 +117,7 @@ void PredictSatellitePeaks::exec() {
   V3D offsets2 = getOffsetVector("ModVector2");
   V3D offsets3 = getOffsetVector("ModVector3");
   int maxOrder = getProperty("MaxOrder");
+  bool crossTerms = getProperty("CrossTerms");
 
   bool includeOrderZero = getProperty("IncludeIntegerHKL");
 
@@ -176,15 +181,21 @@ void PredictSatellitePeaks::exec() {
                                                           true);
   for (auto it = possibleHKLs.begin(); it != possibleHKLs.end(); ++it) {
     V3D hkl = *it;
-    predictOffsets(Peaks, OutPeaks, offsets1, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
-    predictOffsets(Peaks, OutPeaks, offsets2, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
-    predictOffsets(Peaks, OutPeaks, offsets3, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
+    if (crossTerms) {
+      predictOffsetsWithCrossTerms(Peaks, OutPeaks, offsets1, offsets2, offsets3, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+    } else {
+      predictOffsets(Peaks, OutPeaks, 0, offsets1, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+      predictOffsets(Peaks, OutPeaks, 1, offsets2, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+      predictOffsets(Peaks, OutPeaks, 2, offsets3, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+    }
   }
   setProperty("SatellitePeaks", OutPeaks);
 }
@@ -199,6 +210,7 @@ void PredictSatellitePeaks::exec_peaks() {
   V3D offsets2 = getOffsetVector("ModVector2");
   V3D offsets3 = getOffsetVector("ModVector3");
   int maxOrder = getProperty("MaxOrder");
+  bool crossTerms = getProperty("CrossTerms");
 
   bool includePeaksInRange = false;
   bool includeOrderZero = getProperty("IncludeIntegerHKL");
@@ -232,22 +244,28 @@ void PredictSatellitePeaks::exec_peaks() {
   for (auto it = peaks.begin(); it != peaks.end(); ++it) {
     auto peak = *it;
     V3D hkl = peak.getHKL();
-    predictOffsets(Peaks, OutPeaks, offsets1, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
-    predictOffsets(Peaks, OutPeaks, offsets2, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
-    predictOffsets(Peaks, OutPeaks, offsets3, maxOrder, hkl, lambdaFilter,
-                   includePeaksInRange, includeOrderZero, seqNum,
-                   AlreadyDonePeaks);
+    if (crossTerms) {
+      predictOffsetsWithCrossTerms(Peaks, OutPeaks, offsets1, offsets2, offsets3, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+    } else {
+      predictOffsets(Peaks, OutPeaks, 0, offsets1, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+      predictOffsets(Peaks, OutPeaks, 1, offsets2, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+      predictOffsets(Peaks, OutPeaks, 2, offsets3, maxOrder, hkl, lambdaFilter,
+                     includePeaksInRange, includeOrderZero, seqNum,
+                     AlreadyDonePeaks);
+    }
   }
   setProperty("SatellitePeaks", OutPeaks);
 }
 
 void PredictSatellitePeaks::predictOffsets(
     DataObjects::PeaksWorkspace_sptr Peaks,
-    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, V3D offsets,
+    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, int iVector, V3D offsets,
     int &maxOrder, V3D &hkl, HKLFilterWavelength &lambdaFilter,
     bool &includePeaksInRange, bool &includeOrderZero, int &seqNum,
     vector<vector<int>> &AlreadyDonePeaks) {
@@ -297,7 +315,75 @@ void PredictSatellitePeaks::predictOffsets(
     peak->setPeakNumber(seqNum);
     seqNum++;
     peak->setRunNumber(RunNumber);
-    peak->setIntMNP(V3D(order, 0, 0));
+    V3D mnp;
+    mnp[iVector] = order;
+    peak->setIntMNP(mnp);
+    OutPeaks->addPeak(*peak);
+  }
+}
+
+void PredictSatellitePeaks::predictOffsetsWithCrossTerms(
+    DataObjects::PeaksWorkspace_sptr Peaks,
+    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, V3D offsets1,
+    V3D offsets2, V3D offsets3,
+    int &maxOrder, V3D &hkl, HKLFilterWavelength &lambdaFilter,
+    bool &includePeaksInRange, bool &includeOrderZero, int &seqNum,
+    vector<vector<int>> &AlreadyDonePeaks) {
+  if (offsets1 == V3D(0, 0, 0) && offsets2 == V3D(0, 0, 0) && offsets3 == V3D(0, 0, 0))
+    return;
+  const Kernel::DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
+  IPeak &peak1 = Peaks->getPeak(0);
+  Kernel::Matrix<double> goniometer = peak1.getGoniometerMatrix();
+  auto RunNumber = peak1.getRunNumber();
+  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
+  DblMatrix offsetsMat(3, 3);
+  offsetsMat.setRow(0, offsets1);
+  offsetsMat.setRow(1, offsets2);
+  offsetsMat.setRow(2, offsets3);
+  for (int m = -maxOrder; m <= maxOrder; m++)
+    for (int n = -maxOrder; n <= maxOrder; n++)
+      for (int p = -maxOrder; p <= maxOrder; p++) {
+        if (m == 0 && n == 0 && p == 0 && !includeOrderZero)
+          continue; // exclude 0,0,0
+        V3D satelliteHKL(hkl);
+        V3D mnp = V3D(m, n, p);
+        satelliteHKL -= offsetsMat * mnp;
+        if (!lambdaFilter.isAllowed(satelliteHKL) && includePeaksInRange)
+          continue;
+
+
+
+    Kernel::V3D Qs = goniometer * UB * satelliteHKL * 2.0 * M_PI;
+
+    // Check if Q is non-physical
+    if (Qs[2] <= 0)
+      continue;
+
+    auto peak(Peaks->createPeak(Qs, 1));
+
+    peak->setGoniometerMatrix(goniometer);
+
+    if (!peak->findDetector(tracer))
+      continue;
+    vector<int> SavPk{RunNumber, boost::math::iround(1000.0 * satelliteHKL[0]),
+                      boost::math::iround(1000.0 * satelliteHKL[1]),
+                      boost::math::iround(1000.0 * satelliteHKL[2])};
+
+    bool foundPeak =
+        binary_search(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);
+
+    if (!foundPeak) {
+      AlreadyDonePeaks.push_back(SavPk);
+    } else {
+      continue;
+    }
+
+    peak->setHKL(satelliteHKL);
+    peak->setIntHKL(hkl);
+    peak->setPeakNumber(seqNum);
+    seqNum++;
+    peak->setRunNumber(RunNumber);
+    peak->setIntMNP(V3D(m, n, p));
     OutPeaks->addPeak(*peak);
   }
 }