Commit 28c301cb authored by Lynch, Vickie's avatar Lynch, Vickie
Browse files

Refs #22420 function for each offset vector

parent f8aa054e
......@@ -3,6 +3,9 @@
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include "MantidGeometry/Crystal/HKLFilterWavelength.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidAPI/IPeaksWorkspace.h"
namespace Mantid {
namespace Crystal {
......@@ -60,6 +63,7 @@ private:
/// Run the algorithm
void exec() override;
void predictOffsets(DataObjects::PeaksWorkspace_sptr Peaks, boost::shared_ptr<Mantid::API::IPeaksWorkspace>& OutPeaks, std::vector<double> offsets, std::string &label, int &maxOrder, Kernel::V3D &hkl, Kernel::Matrix<double> &goniometer, const Kernel::DblMatrix &UB, Geometry::HKLFilterWavelength &lambdaFilter, bool &includePeaksInRange, bool &includeOrderZero, int &RunNumber, int &seqNum, std::vector<std::vector<int>> &AlreadyDonePeaks);
};
} // namespace Crystal
......
......@@ -289,7 +289,7 @@ namespace Crystal
if (IndexingUtils::ValidIndex(hkl, tolerance))
{
peaks[i].setModulationVector(V3D(0,0,0));
peaks[i].setIntMNP(V3D(0,0,0));
main_indexed++;
h_error = fabs(round(hkl[0])-hkl[0]);
k_error = fabs(round(hkl[1])-hkl[1]);
......@@ -307,7 +307,7 @@ namespace Crystal
hkl1[2] -= order * offsets1[2];
if (IndexingUtils::ValidIndex(hkl1, satetolerance))
{
peaks[i].setModulationVector(V3D(order,0,0));
peaks[i].setIntMNP(V3D(order,0,0));
sate_indexed++;
h_error = fabs(round(hkl1[0])-hkl1[0]);
k_error = fabs(round(hkl1[1])-hkl1[1]);
......@@ -324,7 +324,7 @@ namespace Crystal
hkl1[2] -= order * offsets2[2];
if (IndexingUtils::ValidIndex(hkl1, satetolerance))
{
peaks[i].setModulationVector(V3D(0,order,0));
peaks[i].setIntMNP(V3D(0,order,0));
sate_indexed++;
h_error = fabs(round(hkl1[0])-hkl1[0]);
k_error = fabs(round(hkl1[1])-hkl1[1]);
......@@ -341,7 +341,7 @@ namespace Crystal
hkl1[2] -= order * offsets3[2];
if (IndexingUtils::ValidIndex(hkl1, satetolerance))
{
peaks[i].setModulationVector(V3D(0,0,order));
peaks[i].setIntMNP(V3D(0,0,order));
sate_indexed++;
h_error = fabs(round(hkl1[0])-hkl1[0]);
k_error = fabs(round(hkl1[1])-hkl1[1]);
......
......@@ -393,7 +393,7 @@ DataObjects::Peak LoadIsawPeaks::readPeak(PeaksWorkspace_sptr outWS,
peak.setHKL(h, k, l);
peak.setIntHKL(intHKL);
if (m_isModulatedStructure) {
peak.setModulationVector(mod);
peak.setIntMNP(mod);
}
peak.setIntensity(Inti);
peak.setSigmaIntensity(SigI);
......
......@@ -15,7 +15,6 @@
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidAPI/Run.h"
#include "MantidGeometry/Crystal/BasicHKLFilters.h"
#include "MantidGeometry/Crystal/HKLFilterWavelength.h"
#include "MantidGeometry/Crystal/HKLGenerator.h"
#include <boost/math/special_functions/round.hpp>
......@@ -141,12 +140,6 @@ void PredictSatellitePeaks::exec() {
auto OutPeaks = boost::dynamic_pointer_cast<IPeaksWorkspace>(
WorkspaceFactory::Instance().createPeaks());
OutPeaks->setInstrument(instrument);
OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset1", offsets1,
true);
OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset2", offsets2,
true);
OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset3", offsets3,
true);
V3D hkl;
int peakNum = 0;
......@@ -194,13 +187,48 @@ void PredictSatellitePeaks::exec() {
Progress prog(this, 0.0, 1.0, N);
vector<vector<int>> AlreadyDonePeaks;
bool done = false;
int ErrPos = 1; // Used to determine position in code of a throw
Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
DblMatrix orientedUB = goniometer * UB;
HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
int seqNum = 0;
size_t next = 0;
while (!done) {
std::string offsetName = "Offset1";
predictOffsets(Peaks, OutPeaks, offsets1, offsetName, maxOrder, hkl, goniometer, UB, lambdaFilter, includePeaksInRange, includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
offsetName = "Offset2";
predictOffsets(Peaks, OutPeaks, offsets2, offsetName, maxOrder, hkl, goniometer, UB, lambdaFilter, includePeaksInRange, includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
offsetName = "Offset3";
predictOffsets(Peaks, OutPeaks, offsets3, offsetName, maxOrder, hkl, goniometer, UB, lambdaFilter, includePeaksInRange, includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
if (includePeaksInRange) {
next++;
if (next == possibleHKLs.size())
break;
hkl = possibleHKLs[next];
} 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();
goniometer = peak1.getGoniometerMatrix();
RunNumber = peak1.getRunNumber();
}
}
prog.report();
}
setProperty("SatellitePeaks", OutPeaks);
}
void PredictSatellitePeaks::predictOffsets(DataObjects::PeaksWorkspace_sptr Peaks, boost::shared_ptr<Mantid::API::IPeaksWorkspace>& OutPeaks, std::vector<double> offsets, std::string &label, int &maxOrder, V3D &hkl, Kernel::Matrix<double> &goniometer, const Kernel::DblMatrix &UB, HKLFilterWavelength &lambdaFilter, bool &includePeaksInRange, bool &includeOrderZero, int &RunNumber, int &seqNum, vector<vector<int>> &AlreadyDonePeaks) {
OutPeaks->mutableRun().addProperty<std::vector<double>>(label, offsets,
true);
int ErrPos = 1; // Used to determine position in code of a throw
Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
for (int order = -maxOrder; order <= maxOrder; order++) {
if (order == 0 && !includeOrderZero)
continue; // exclude order 0
......@@ -208,9 +236,9 @@ void PredictSatellitePeaks::exec() {
V3D hkl1(hkl);
ErrPos = 0;
hkl1[0] += order * offsets1[0];
hkl1[1] += order * offsets1[1];
hkl1[2] += order * offsets1[2];
hkl1[0] += order * offsets[0];
hkl1[1] += order * offsets[1];
hkl1[2] += order * offsets[2];
if (!lambdaFilter.isAllowed(hkl1) && includePeaksInRange)
continue;
......@@ -245,7 +273,7 @@ void PredictSatellitePeaks::exec() {
peak->setPeakNumber(seqNum);
seqNum++;
peak->setRunNumber(RunNumber);
peak->setModulationVector(V3D(order, 0, 0));
peak->setIntMNP(V3D(order, 0, 0));
OutPeaks->addPeak(*peak);
}
} catch (std::runtime_error &) {
......@@ -253,127 +281,8 @@ void PredictSatellitePeaks::exec() {
throw std::invalid_argument("Invalid data at this point");
}
}
for (int order = -maxOrder; order <= maxOrder; 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 = goniometer * UB * hkl1 * 2.0 * M_PI;
if (Qs[2] <= 0)
continue;
ErrPos = 1;
boost::shared_ptr<IPeak> peak(Peaks->createPeak(Qs, 1));
peak->setGoniometerMatrix(goniometer);
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->setPeakNumber(seqNum);
seqNum++;
peak->setRunNumber(RunNumber);
peak->setModulationVector(V3D(0, order, 0));
OutPeaks->addPeak(*peak);
}
} catch (std::runtime_error &) {
if (ErrPos != 1) // setQLabFrame in createPeak throws exception
throw std::invalid_argument("Invalid data at this point");
}
}
for (int order = -maxOrder; order <= maxOrder; 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 = goniometer * UB * hkl1 * 2.0 * M_PI;
if (Qs[2] <= 0)
continue;
ErrPos = 1;
boost::shared_ptr<IPeak> peak(Peaks->createPeak(Qs, 1));
peak->setGoniometerMatrix(goniometer);
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->setPeakNumber(seqNum);
seqNum++;
peak->setRunNumber(RunNumber);
peak->setModulationVector(V3D(0, 0, order));
OutPeaks->addPeak(*peak);
}
} catch (std::runtime_error &) {
if (ErrPos != 1) // setQLabFrame in createPeak throws exception
throw std::invalid_argument("Invalid data at this point");
}
}
if (includePeaksInRange) {
next++;
if (next == possibleHKLs.size())
break;
hkl = possibleHKLs[next];
} 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();
goniometer = peak1.getGoniometerMatrix();
RunNumber = peak1.getRunNumber();
}
}
prog.report();
}
setProperty("SatellitePeaks", OutPeaks);
}
} // namespace Crystal
} // namespace Mantid
......@@ -94,7 +94,7 @@ void SaveIsawPeaks::exec() {
runMap_t runMap;
for (size_t i = 0; i < peaks.size(); ++i) {
Peak &p = peaks[i];
if (p.getModulationVector() != V3D(0, 0, 0))
if (p.getIntMNP() != V3D(0, 0, 0))
m_isModulatedStructure = true;
int run = p.getRunNumber();
int bank = 0;
......@@ -353,7 +353,7 @@ void SaveIsawPeaks::exec() {
// HKL's are flipped by -1 because of the internal Q convention
// unless Crystallography convention
if (m_isModulatedStructure) {
V3D mod = p.getModulationVector();
V3D mod = p.getIntMNP();
auto intHKL = p.getIntHKL();
out << std::setw(5) << Utils::round(qSign * intHKL.X())
<< std::setw(5) << Utils::round(qSign * intHKL.Y())
......
......@@ -154,8 +154,8 @@ public:
void setCol(int m_col);
void setPeakNumber(int m_peakNumber) override;
int getPeakNumber() const override;
void setModulationVector(const Mantid::Kernel::V3D m_modStru) override;
Mantid::Kernel::V3D getModulationVector() const override;
void setIntMNP(const Mantid::Kernel::V3D m_modStru) override;
Mantid::Kernel::V3D getIntMNP() const override;
virtual Mantid::Kernel::V3D getDetPos() const override;
double getL1() const override;
......
......@@ -229,7 +229,7 @@ Peak::Peak(const Geometry::IPeak &ipeak)
m_runNumber(ipeak.getRunNumber()),
m_monitorCount(ipeak.getMonitorCount()), m_row(ipeak.getRow()),
m_col(ipeak.getCol()), m_orig_H(0.), m_orig_K(0.), m_orig_L(0.),
m_peakNumber(ipeak.getPeakNumber()), m_modStru(ipeak.getModulationVector()),
m_peakNumber(ipeak.getPeakNumber()), m_modStru(ipeak.getIntMNP()),
m_intH(int(ipeak.getIntHKL().X())), m_intK(int(ipeak.getIntHKL().Y())),
m_intL(int(ipeak.getIntHKL().Z())),
m_peakShape(boost::make_shared<NoShape>()) {
......@@ -924,7 +924,7 @@ int Peak::getPeakNumber() const { return m_peakNumber; }
// -------------------------------------------------------------------------------------
/**Returns the unique peak number
* Returns -1 if it could not find it. */
V3D Peak::getModulationVector() const { return m_modStru; }
V3D Peak::getIntMNP() const { return m_modStru; }
// -------------------------------------------------------------------------------------
/** For RectangularDetectors only, sets the row (y) of the pixel of the
......@@ -948,7 +948,7 @@ void Peak::setPeakNumber(int m_peakNumber) {
// -------------------------------------------------------------------------------------
/** Sets the modulated peak structure number
* @param m_modStru :: modulated peak structure value */
void Peak::setModulationVector(V3D m_modStru) { this->m_modStru = m_modStru; }
void Peak::setIntMNP(V3D m_modStru) { this->m_modStru = m_modStru; }
// -------------------------------------------------------------------------------------
/** Return the detector position vector */
......
......@@ -84,8 +84,8 @@ public:
virtual int getPeakNumber() const = 0;
virtual void setPeakNumber(int m_PeakNumber) = 0;
virtual Mantid::Kernel::V3D getModulationVector() const = 0;
virtual void setModulationVector(const Mantid::Kernel::V3D m_modStru) = 0;
virtual Mantid::Kernel::V3D getIntMNP() const = 0;
virtual void setIntMNP(const Mantid::Kernel::V3D m_modStru) = 0;
virtual Mantid::Kernel::Matrix<double> getGoniometerMatrix() const = 0;
virtual void setGoniometerMatrix(
......
......@@ -67,10 +67,10 @@ public:
MOCK_CONST_METHOD0(getInstrument, Geometry::Instrument_const_sptr());
MOCK_CONST_METHOD0(getRunNumber, int());
MOCK_CONST_METHOD0(getPeakNumber, int());
MOCK_CONST_METHOD0(getModulationVector, Mantid::Kernel::V3D());
MOCK_CONST_METHOD0(getIntMNP, Mantid::Kernel::V3D());
MOCK_METHOD1(setRunNumber, void(int m_RunNumber));
MOCK_METHOD1(setPeakNumber, void(int m_PeakNumber));
MOCK_METHOD1(setModulationVector, void(const Mantid::Kernel::V3D m_modStru));
MOCK_METHOD1(setIntMNP, void(const Mantid::Kernel::V3D m_modStru));
MOCK_CONST_METHOD0(getMonitorCount, double());
MOCK_METHOD1(setMonitorCount, void(double m_MonitorCount));
MOCK_CONST_METHOD0(getH, double());
......
......@@ -61,7 +61,7 @@ void export_IPeak() {
"cache values related to it.")
.def("getRunNumber", &IPeak::getRunNumber, arg("self"),
"Return the run number this peak was measured at")
.def("getModulationVector", &IPeak::getModulationVector, arg("self"),
.def("getIntMNP", &IPeak::getIntMNP, arg("self"),
"Return the modulated scructure for this peak")
.def("getPeakNumber", &IPeak::getPeakNumber, arg("self"),
"Return the peak number for this peak")
......@@ -70,7 +70,7 @@ void export_IPeak() {
.def("setRunNumber", &IPeak::setRunNumber,
(arg("self"), arg("run_number")),
"Set the run number that measured this peak")
.def("setModulationVector", &IPeak::setModulationVector,
.def("setIntMNP", &IPeak::setIntMNP,
(arg("self"), arg("modulated_structure")),
"Set the modulated structure for this peak")
.def("setPeakNumber", &IPeak::setPeakNumber,
......
......@@ -183,8 +183,8 @@ class IPeakTest(unittest.TestCase):
def test_set_modulation_vector(self):
testVector = V3D(0.5,0,0.2)
self._peak.setModulationVector(testVector)
self.assertEqual(self._peak.getModulationVector(), testVector)
self._peak.setIntMNP(testVector)
self.assertEqual(self._peak.getIntMNP(), testVector)
if __name__ == '__main__':
unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment