Refs #22420 add order columns to PredictSatellitePeaks

Framework/Crystal/CMakeLists.txt

@@ -51,6 +51,7 @@ set ( SRC_FILES
 	src/PeaksOnSurface.cpp
 	src/PredictFractionalPeaks.cpp
 	src/PredictPeaks.cpp
+	src/PredictSatellitePeaks.cpp
 	src/SCDCalibratePanels.cpp
 	src/SCDPanelErrors.cpp
 	src/SaveHKL.cpp
@@ -127,6 +128,7 @@ set ( INC_FILES
 	inc/MantidCrystal/PeaksOnSurface.h
 	inc/MantidCrystal/PredictFractionalPeaks.h
 	inc/MantidCrystal/PredictPeaks.h
+	inc/MantidCrystal/PredictSatellitePeaks.h
 	inc/MantidCrystal/SCDCalibratePanels.h
 	inc/MantidCrystal/SCDPanelErrors.h
 	inc/MantidCrystal/SaveHKL.h
@@ -196,6 +198,7 @@ set ( TEST_FILES
 	PeaksOnSurfaceTest.h
 	PredictFractionalPeaksTest.h
 	PredictPeaksTest.h
+	PredictSatellitePeaksTest.h
 	SCDCalibratePanelsTest.h
 	SaveHKLTest.h
 	SaveIsawPeaksTest.h

Framework/Crystal/inc/MantidCrystal/PredictSatellitePeaks.h

+#ifndef MANTID_CRYSTAL_PREDICTSATELLITEPEAKS_H_
+#define MANTID_CRYSTAL_PREDICTSATELLITEPEAKS_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid {
+namespace Crystal {
+/** CreatSatellitePeaks : Algorithm to create a PeaksWorkspace with peaks
+   corresponding
+    to fractional h,k,and l values.
+
+    @author Ruth Mikkelson
+    @date   2012-12-05
+
+    Copyright © 2012 ISIS Rutherford Appleton Laboratory &
+                     NScD Oak Ridge National Laboratory
+
+    This file is part of Mantid.
+
+    Mantid is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.
+
+    Mantid is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+    File change history is stored at:
+    <https://github.com/mantidproject/mantid>
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+  */
+class DLLExport PredictSatellitePeaks : public API::Algorithm {
+public:
+  /// Algorithm's name for identification
+  const std::string name() const override { return "PredictSatellitePeaks"; };
+  /// Summary of algorithms purpose
+  const std::string summary() const override {
+    return "The offsets can be from hkl values in a range of hkl values or "
+           "from peaks in the input PeaksWorkspace";
+  }
+
+  /// Algorithm's version for identification
+  int version() const override { return 1; };
+  const std::vector<std::string> seeAlso() const override {
+    return {"PredictPeaks"};
+  }
+
+  /// Algorithm's category for identification
+  const std::string category() const override { return "Crystal\\Peaks"; }
+
+private:
+  /// Initialise the properties
+  void init() override;
+
+  /// Run the algorithm
+  void exec() override;
+};
+
+} // namespace Crystal
+} // namespace Mantid
+
+#endif /* MANTID_CRYSTAL_PREDICTSATELLITEPEAKS */

Framework/Crystal/src/PredictSatellitePeaks.cpp

+/*
+* PredictSatellitePeaks.cpp
+*
+*  Created on: Dec 5, 2012
+*      Author: ruth
+*/
+#include "MantidCrystal/PredictSatellitePeaks.h"
+#include "MantidAPI/Sample.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+#include "MantidGeometry/Crystal/OrientedLattice.h"
+#include "MantidGeometry/Objects/InstrumentRayTracer.h"
+#include "MantidKernel/ArrayProperty.h"
+#include "MantidKernel/ArrayLengthValidator.h"
+#include "MantidKernel/EnabledWhenProperty.h"
+
+#include <boost/math/special_functions/round.hpp>
+
+namespace Mantid {
+using namespace Mantid::DataObjects;
+using namespace Mantid::API;
+using namespace std;
+using namespace Mantid::Geometry;
+using namespace Mantid::Kernel;
+namespace Crystal {
+
+DECLARE_ALGORITHM(PredictSatellitePeaks)
+
+/// Initialise the properties
+void PredictSatellitePeaks::init() {
+  declareProperty(
+      make_unique<WorkspaceProperty<PeaksWorkspace>>("Peaks", "",
+                                                     Direction::Input),
+      "Workspace of Peaks with orientation matrix that indexed the peaks and "
+      "instrument loaded");
+
+  declareProperty(
+      make_unique<WorkspaceProperty<PeaksWorkspace>>("SatellitePeaks", "",
+                                                     Direction::Output),
+      "Workspace of Peaks with peaks with fractional h,k, and/or l values");
+  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
+                      string("OffsetVector1"), "0.0,0.0,0.0,0"),
+                  "Offsets for h, k, l directions and order");
+  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
+                      string("OffsetVector2"), "0.0,0.0,0.0,0"),
+                  "Offsets for h, k, l directions and order");
+  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
+                      string("OffsetVector3"), "0.0,0.0,0.0,0"),
+                  "Offsets for h, k, l directions and order");
+
+  declareProperty("IncludeAllPeaksInRange", false,
+                  "If false only offsets from peaks from Peaks are used");
+
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Hmin", -8.0, Direction::Input),
+      "Minimum H value to use");
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Hmax", 8.0, Direction::Input),
+      "Maximum H value to use");
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Kmin", -8.0, Direction::Input),
+      "Minimum K value to use");
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Kmax", 8.0, Direction::Input),
+      "Maximum K value to use");
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Lmin", -8.0, Direction::Input),
+      "Minimum L value to use");
+  declareProperty(
+      make_unique<PropertyWithValue<double>>("Lmax", 8.0, Direction::Input),
+      "Maximum L value to use");
+
+  setPropertySettings(
+      "Hmin", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+
+  setPropertySettings(
+      "Hmax", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+  setPropertySettings(
+      "Kmin", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+
+  setPropertySettings(
+      "Kmax", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+  setPropertySettings(
+      "Lmin", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+
+  setPropertySettings(
+      "Lmax", Kernel::make_unique<Kernel::EnabledWhenProperty>(
+                  string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
+}
+
+/// Run the algorithm
+void PredictSatellitePeaks::exec() {
+  PeaksWorkspace_sptr Peaks = getProperty("Peaks");
+  if (!Peaks)
+    throw std::invalid_argument(
+        "Input workspace is not a PeaksWorkspace. Type=" + Peaks->id());
+
+  vector<double> offsets1 = getProperty("OffsetVector1");
+  vector<double> offsets2 = getProperty("OffsetVector2");
+  vector<double> offsets3 = getProperty("OffsetVector3");
+  if (offsets1.empty()) {
+    offsets1.push_back(0.0);
+    offsets1.push_back(0.0);
+    offsets1.push_back(0.0);
+    offsets1.push_back(0.0);
+  }
+  if (offsets2.empty()) {
+    offsets2.push_back(0.0);
+    offsets2.push_back(0.0);
+    offsets2.push_back(0.0);
+    offsets2.push_back(0.0);
+  }
+  if (offsets3.empty()) {
+    offsets3.push_back(0.0);
+    offsets3.push_back(0.0);
+    offsets3.push_back(0.0);
+    offsets3.push_back(0.0);
+  }
+
+  bool includePeaksInRange = getProperty("IncludeAllPeaksInRange");
+
+  if (Peaks->getNumberPeaks() <= 0) {
+    g_log.error() << "There are No peaks in the input PeaksWorkspace\n";
+    return;
+  }
+
+  API::Sample samp = Peaks->sample();
+
+  Geometry::OrientedLattice &ol = samp.getOrientedLattice();
+
+  Geometry::Instrument_const_sptr Instr = Peaks->getInstrument();
+
+  auto OutPeaks = boost::dynamic_pointer_cast<IPeaksWorkspace>(
+      WorkspaceFactory::Instance().createPeaks());
+  OutPeaks->setInstrument(Instr);
+
+  V3D hkl;
+  int peakNum = 0;
+  const auto NPeaks = Peaks->getNumberPeaks();
+  Kernel::Matrix<double> Gon;
+  Gon.identityMatrix();
+
+  const double Hmin = getProperty("Hmin");
+  const double Hmax = getProperty("Hmax");
+  const double Kmin = getProperty("Kmin");
+  const double Kmax = getProperty("Kmax");
+  const double Lmin = getProperty("Lmin");
+  const double Lmax = getProperty("Lmax");
+
+  int N = NPeaks;
+  if (includePeaksInRange) {
+    N = boost::math::iround((Hmax - Hmin + 1) * (Kmax - Kmin + 1) *
+                            (Lmax - Lmin + 1));
+    N = max<int>(100, N);
+  }
+  IPeak &peak0 = Peaks->getPeak(0);
+  auto RunNumber = peak0.getRunNumber();
+  Gon = peak0.getGoniometerMatrix();
+  Progress prog(this, 0.0, 1.0, N);
+  if (includePeaksInRange) {
+
+    hkl[0] = Hmin;
+    hkl[1] = Kmin;
+    hkl[2] = Lmin;
+  } else {
+    hkl[0] = peak0.getH();
+    hkl[1] = peak0.getK();
+    hkl[2] = peak0.getL();
+  }
+
+  const Kernel::DblMatrix &UB = ol.getUB();
+  vector<vector<int>> AlreadyDonePeaks;
+  bool done = false;
+  int ErrPos = 1; // Used to determine position in code of a throw
+  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
+  while (!done) {
+    for (int order = -static_cast<int>(offsets1[3]); order <= static_cast<int>(offsets1[3]); order++) {
+          try {
+            V3D hkl1(hkl);
+            ErrPos = 0;
+
+            hkl1[0] += order * offsets1[0];
+            hkl1[1] += order * offsets1[1];
+            hkl1[2] += order * offsets1[2];
+
+            Kernel::V3D Qs = UB * hkl1;
+            Qs *= 2.0;
+            Qs *= M_PI;
+            Qs = Gon * Qs;
+            if (Qs[2] <= 0)
+              continue;
+
+            ErrPos = 1;
+
+            boost::shared_ptr<IPeak> peak(Peaks->createPeak(Qs, 1));
+
+            peak->setGoniometerMatrix(Gon);
+
+            if (Qs[2] > 0 && peak->findDetector(tracer)) {
+              ErrPos = 2;
+              vector<int> SavPk{RunNumber,
+                                boost::math::iround(1000.0 * hkl1[0]),
+                                boost::math::iround(1000.0 * hkl1[1]),
+                                boost::math::iround(1000.0 * hkl1[2])};
+
+              // TODO keep list sorted so searching is faster?
+              auto it =
+                  find(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);
+
+              if (it == AlreadyDonePeaks.end())
+                AlreadyDonePeaks.push_back(SavPk);
+              else
+                continue;
+
+              peak->setHKL(hkl1);
+              peak->setRunNumber(RunNumber);
+              peak->setModStru(V3D(order,0,0));
+              OutPeaks->addPeak(*peak);
+            }
+          } catch (...) {
+            if (ErrPos != 1) // setQLabFrame in createPeak throws exception
+              throw std::invalid_argument("Invalid data at this point");
+    }
+    }
+    for (int order = -static_cast<int>(offsets2[3]); order <= static_cast<int>(offsets2[3]); order++) {
+        if (order == 0) continue; // already added with 1st vector
+          try {
+            V3D hkl1(hkl);
+            ErrPos = 0;
+
+            hkl1[0] += order * offsets2[0];
+            hkl1[1] += order * offsets2[1];
+            hkl1[2] += order * offsets2[2];
+
+            Kernel::V3D Qs = UB * hkl1;
+            Qs *= 2.0;
+            Qs *= M_PI;
+            Qs = Gon * Qs;
+            if (Qs[2] <= 0)
+              continue;
+
+            ErrPos = 1;
+
+            boost::shared_ptr<IPeak> peak(Peaks->createPeak(Qs, 1));
+
+            peak->setGoniometerMatrix(Gon);
+
+            if (Qs[2] > 0 && peak->findDetector(tracer)) {
+              ErrPos = 2;
+              vector<int> SavPk{RunNumber,
+                                boost::math::iround(1000.0 * hkl1[0]),
+                                boost::math::iround(1000.0 * hkl1[1]),
+                                boost::math::iround(1000.0 * hkl1[2])};
+
+              // TODO keep list sorted so searching is faster?
+              auto it =
+                  find(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);
+
+              if (it == AlreadyDonePeaks.end())
+                AlreadyDonePeaks.push_back(SavPk);
+              else
+                continue;
+
+              peak->setHKL(hkl1);
+              peak->setRunNumber(RunNumber);
+              peak->setModStru(V3D(0,order,0));
+              OutPeaks->addPeak(*peak);
+            }
+          } catch (...) {
+            if (ErrPos != 1) // setQLabFrame in createPeak throws exception
+              throw std::invalid_argument("Invalid data at this point");
+    }
+    }
+    for (int order = -static_cast<int>(offsets3[3]); order <= static_cast<int>(offsets3[3]); order++) {
+        if (order == 0) continue; // already added with 1st vector
+          try {
+            V3D hkl1(hkl);
+            ErrPos = 0;
+
+            hkl1[0] += order * offsets3[0];
+            hkl1[1] += order * offsets3[1];
+            hkl1[2] += order * offsets3[2];
+
+            Kernel::V3D Qs = UB * hkl1;
+            Qs *= 2.0;
+            Qs *= M_PI;
+            Qs = Gon * Qs;
+            if (Qs[2] <= 0)
+              continue;
+
+            ErrPos = 1;
+
+            boost::shared_ptr<IPeak> peak(Peaks->createPeak(Qs, 1));
+
+            peak->setGoniometerMatrix(Gon);
+
+            if (Qs[2] > 0 && peak->findDetector(tracer)) {
+              ErrPos = 2;
+              vector<int> SavPk{RunNumber,
+                                boost::math::iround(1000.0 * hkl1[0]),
+                                boost::math::iround(1000.0 * hkl1[1]),
+                                boost::math::iround(1000.0 * hkl1[2])};
+
+              // TODO keep list sorted so searching is faster?
+              auto it =
+                  find(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);
+
+              if (it == AlreadyDonePeaks.end())
+                AlreadyDonePeaks.push_back(SavPk);
+              else
+                continue;
+
+              peak->setHKL(hkl1);
+              peak->setRunNumber(RunNumber);
+              peak->setModStru(V3D(0,0,order));
+              OutPeaks->addPeak(*peak);
+            }
+          } catch (...) {
+            if (ErrPos != 1) // setQLabFrame in createPeak throws exception
+              throw std::invalid_argument("Invalid data at this point");
+    }
+    }
+    if (includePeaksInRange) {
+      hkl[0]++;
+      if (hkl[0] > Hmax) {
+        hkl[0] = Hmin;
+        hkl[1]++;
+        if (hkl[1] > Kmax) {
+
+          hkl[1] = Kmin;
+          hkl[2]++;
+          if (hkl[2] > Lmax)
+            done = true;
+        }
+      }
+    } else {
+      peakNum++;
+      if (peakNum >= NPeaks)
+        done = true;
+      else { // peak0= Peaks->getPeak(peakNum);
+        IPeak &peak1 = Peaks->getPeak(peakNum);
+        //??? could not assign to peak0 above. Did not work
+        // the peak that peak0 was associated with did NOT change
+        hkl[0] = peak1.getH();
+        hkl[1] = peak1.getK();
+        hkl[2] = peak1.getL();
+        Gon = peak1.getGoniometerMatrix();
+        RunNumber = peak1.getRunNumber();
+      }
+    }
+    prog.report();
+  }
+
+  setProperty("SatellitePeaks", OutPeaks);
+}
+
+} // namespace Crystal
+} // namespace Mantid

Framework/Crystal/src/SaveIsawPeaks.cpp

@@ -53,6 +53,9 @@ void SaveIsawPeaks::exec() {
   std::string header = "2   SEQN    H    K    L     COL      ROW     CHAN      "
                        "  L2   2_THETA        AZ         WL         D      IPK "
                        "      INTI    SIGI  RFLG";
+  header = "2   SEQN    H    K    L     ORD1 ORD2 ORD3  COL      ROW     CHAN      "
+                       "  L2   2_THETA        AZ         WL         D      IPK "
+                       "      INTI    SIGI  RFLG";
 
   std::string filename = getPropertyValue("Filename");
   PeaksWorkspace_sptr ws = getProperty("InputWorkspace");
@@ -155,7 +158,7 @@ void SaveIsawPeaks::exec() {
     // saving.
     Types::Core::DateAndTime expDate = inst->getValidFromDate() + 1.0;
     out << expDate.toISO8601String();
-    if (m_ModStru) out << "MOD";
+    if (m_ModStru) out << " MOD";
     out << '\n';
 
     out << "6         L1    T0_SHIFT\n";
@@ -320,14 +323,18 @@ void SaveIsawPeaks::exec() {
 
           // HKL's are flipped by -1 because of the internal Q convention
           // unless Crystallography convention
+          if (m_ModStru) {
+          V3D mod = p.getModStru();
+          out << std::setw(5) << Utils::round(qSign * (p.getH()-mod[0]*0.5)) << std::setw(5)
+              << Utils::round(qSign * (p.getK()-mod[1])) << std::setw(5)
+              << Utils::round(qSign * (p.getL()-mod[2]));
+             out << std::setw(5) <<  Utils::round(qSign * mod[0]) << std::setw(5)
+              << Utils::round(qSign * mod[1]) << std::setw(5) << Utils::round(qSign * mod[2]);
+          }
+          else {
           out << std::setw(5) << Utils::round(qSign * p.getH()) << std::setw(5)
               << Utils::round(qSign * p.getK()) << std::setw(5)
               << Utils::round(qSign * p.getL());
-
-          if (m_ModStru) {
-          V3D mod = p.getModStru();
-             out << std::setw(5) <<  mod[0] << std::setw(5)
-              << mod[1] << std::setw(5) << mod[2];
           }
 
           // Row/column

Framework/Crystal/test/PredictSatellitePeaksTest.h

+/*
+ * PredictSatellitePeaksTest.h
+ *
+ *  Created on: Dec 28, 2012
+ *      Author: ruth
+ */
+
+#ifndef PREDICTFRACTIONALPEAKSTEST_H_
+#define PREDICTFRACTIONALPEAKSTEST_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidKernel/Timer.h"
+#include "MantidKernel/System.h"
+
+#include "MantidDataHandling/LoadNexusProcessed.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+#include "MantidCrystal/LoadIsawUB.h"
+#include "MantidCrystal/IndexPeaks.h"
+#include "MantidCrystal/PredictSatellitePeaks.h"
+#include "MantidAPI/FrameworkManager.h"
+
+using namespace Mantid::Crystal;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+using namespace Mantid::DataHandling;
+using namespace Mantid::Kernel;
+
+class PredictSatellitePeaksTest : public CxxTest::TestSuite {
+
+public:
+  void test_Init() {
+    PredictSatellitePeaks alg;
+    TS_ASSERT_THROWS_NOTHING(alg.initialize());
+    TS_ASSERT(alg.isInitialized());
+  }
+
+  void test_exec() {
+    LoadNexusProcessed loader;
+    TS_ASSERT_THROWS_NOTHING(loader.initialize());
+    TS_ASSERT(loader.isInitialized());
+
+    std::string WSName("peaks");
+    loader.setPropertyValue("Filename", "TOPAZ_3007.peaks.nxs");
+    loader.setPropertyValue("OutputWorkspace", WSName);
+    TS_ASSERT(loader.execute());
+    TS_ASSERT(loader.isExecuted());
+
+    LoadIsawUB UBloader;
+    TS_ASSERT_THROWS_NOTHING(UBloader.initialize());
+    TS_ASSERT(UBloader.isInitialized());
+
+    UBloader.setPropertyValue("InputWorkspace", WSName);
+
+    UBloader.setProperty("FileName", "TOPAZ_3007.mat");
+    TS_ASSERT(UBloader.execute());
+    TS_ASSERT(UBloader.isExecuted());
+
+    IndexPeaks indexer;
+    TS_ASSERT_THROWS_NOTHING(indexer.initialize());
+    TS_ASSERT(indexer.isInitialized());
+
+    indexer.setPropertyValue("PeaksWorkspace", WSName);
+
+    TS_ASSERT(indexer.execute());
+    TS_ASSERT(indexer.isExecuted());
+
+    PredictSatellitePeaks alg;
+    TS_ASSERT_THROWS_NOTHING(alg.initialize());
+    TS_ASSERT(alg.isInitialized());
+
+    alg.setProperty("Peaks", WSName);
+
+    alg.setProperty("FracPeaks", std::string("FracPeaks"));
+    alg.setProperty("Offset1", "0.5,0,.2,1");
+    TS_ASSERT(alg.execute());
+    TS_ASSERT(alg.isExecuted());
+    alg.setPropertyValue("FracPeaks", "FracPeaks");
+    PeaksWorkspace_sptr FracPeaks = alg.getProperty("FracPeaks");
+
+    TS_ASSERT_EQUALS(FracPeaks->getNumberPeaks(), 117);
+
+    Peak &peakx = dynamic_cast<Peak &>(FracPeaks->getPeak(0));
+    Peak peak = peakx;
+    TS_ASSERT_DELTA(peak.getH(), -5.5, .0001);
+    TS_ASSERT_DELTA(peak.getK(), 7.0, .0001);
+    TS_ASSERT_DELTA(peak.getL(), -3.8, .0001);
+
+    peakx = dynamic_cast<Peak &>(FracPeaks->getPeak(3));
+    peak = peakx;
+    TS_ASSERT_DELTA(peak.getH(), -5.5, .0001);
+    TS_ASSERT_DELTA(peak.getK(), 3.0, .0001);
+    TS_ASSERT_DELTA(peak.getL(), -2.8, .0001);
+
+    peakx = dynamic_cast<Peak &>(FracPeaks->getPeak(6));
+    peak = peakx;
+    TS_ASSERT_DELTA(peak.getH(), -6.5, .0001);