Unverified Commit 397e0ca9 authored by Whitfield, Ross's avatar Whitfield, Ross Committed by GitHub
Browse files

Merge pull request #31227 from rosswhitfield/demand_satellite_peaks

Update DEMAND workflow to support satellite peaks
parents 85255fd1 b9629131
......@@ -11,7 +11,6 @@
#include "MantidAPI/ILatticeFunction.h"
#include "MantidCrystal/DllConfig.h"
#include "MantidDataObjects/OffsetsWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/UnitCell.h"
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
......
......@@ -6,6 +6,7 @@
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCrystal/CalculateUMatrix.h"
#include "MantidAPI/Sample.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/BoundedValidator.h"
......@@ -24,7 +25,7 @@ using Mantid::Geometry::OrientedLattice;
/** Initialize the algorithm's properties.
*/
void CalculateUMatrix::init() {
this->declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
"An input workspace.");
std::shared_ptr<BoundedValidator<double>> mustBePositive = std::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
......@@ -53,7 +54,7 @@ void CalculateUMatrix::exec() {
auto lattice = std::make_unique<OrientedLattice>(a, b, c, alpha, beta, gamma);
Matrix<double> B = lattice->getB();
PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
IPeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
if (!ws)
throw std::runtime_error("Problems reading the peaks workspace");
......@@ -62,7 +63,7 @@ void CalculateUMatrix::exec() {
V3D old(0, 0, 0);
Matrix<double> Hi(4, 4), Si(4, 4), HS(4, 4), zero(4, 4);
for (int i = 0; i < ws->getNumberPeaks(); i++) {
Peak p = ws->getPeaks()[i];
Mantid::Geometry::IPeak &p = ws->getPeak(i);
double H = p.getH();
double K = p.getK();
double L = p.getL();
......
......@@ -10,7 +10,10 @@
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidCrystal/GSLFunctions.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/EdgePixel.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
......@@ -37,7 +40,7 @@ using namespace DataObjects;
*/
void OptimizeLatticeForCellType::init() {
declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
"An input PeaksWorkspace with an instrument.");
std::vector<std::string> cellTypes;
cellTypes.emplace_back(ReducedCell::CUBIC());
......@@ -73,20 +76,21 @@ void OptimizeLatticeForCellType::exec() {
double tolerance = this->getProperty("Tolerance");
int edge = this->getProperty("EdgePixels");
std::string cell_type = getProperty("CellType");
DataObjects::PeaksWorkspace_sptr ws = getProperty("PeaksWorkspace");
Geometry::Instrument_const_sptr inst = ws->getInstrument();
IPeaksWorkspace_sptr ws = getProperty("PeaksWorkspace");
std::vector<int> badPeaks;
std::vector<DataObjects::PeaksWorkspace_sptr> runWS;
if (edge > 0) {
for (int i = int(ws->getNumberPeaks()) - 1; i >= 0; --i) {
const std::vector<Peak> &peaks = ws->getPeaks();
if (edgePixel(inst, peaks[i].getBankName(), peaks[i].getCol(), peaks[i].getRow(), edge)) {
badPeaks.emplace_back(i);
std::vector<IPeaksWorkspace_sptr> runWS;
if (edge > 0)
if (auto pw = std::dynamic_pointer_cast<PeaksWorkspace>(ws)) {
Geometry::Instrument_const_sptr inst = ws->getInstrument();
for (int i = int(pw->getNumberPeaks()) - 1; i >= 0; --i) {
const std::vector<Peak> &peaks = pw->getPeaks();
if (edgePixel(inst, peaks[i].getBankName(), peaks[i].getCol(), peaks[i].getRow(), edge)) {
badPeaks.emplace_back(i);
}
}
pw->removePeaks(std::move(badPeaks));
}
ws->removePeaks(std::move(badPeaks));
}
runWS.emplace_back(ws);
if (perRun) {
......@@ -94,14 +98,13 @@ void OptimizeLatticeForCellType::exec() {
// Sort by run number
criteria.emplace_back("runnumber", true);
ws->sort(criteria);
const std::vector<Peak> &peaks_all = ws->getPeaks();
int run = 0;
int count = 0;
for (const auto &peak : peaks_all) {
for (int i = 0; i < ws->getNumberPeaks(); i++) {
IPeak &peak = ws->getPeak(i);
if (peak.getRunNumber() != run) {
count++; // first entry in runWS is input workspace
auto cloneWS = std::make_shared<PeaksWorkspace>();
cloneWS->setInstrument(inst);
auto cloneWS = std::dynamic_pointer_cast<IPeaksWorkspace>(WorkspaceFactory::Instance().createPeaks(ws->id()));
cloneWS->copyExperimentInfoFrom(ws.get());
runWS.emplace_back(cloneWS);
runWS[count]->addPeak(peak);
......@@ -114,7 +117,7 @@ void OptimizeLatticeForCellType::exec() {
}
// finally do the optimization
for (auto &i_run : runWS) {
DataObjects::PeaksWorkspace_sptr peakWS(i_run->clone());
IPeaksWorkspace_sptr peakWS(i_run->clone());
AnalysisDataService::Instance().addOrReplace("_peaks", peakWS);
const DblMatrix UB = peakWS->sample().getOrientedLattice().getUB();
auto ol = peakWS->sample().getOrientedLattice();
......
......@@ -21,8 +21,6 @@
#include "MantidGeometry/Crystal/HKLGenerator.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include <boost/math/special_functions/round.hpp>
......@@ -53,19 +51,10 @@ void PredictSatellitePeaks::init() {
declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("SatellitePeaks", "", Direction::Output),
"Workspace of Peaks with peaks with fractional h,k, and/or l values");
declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(string("ModVector1"), "0.0,0.0,0.0"),
"Offsets for h, k, l directions ");
declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(string("ModVector2"), "0.0,0.0,0.0"),
"Offsets for h, k, l directions ");
declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(string("ModVector3"), "0.0,0.0,0.0"),
"Offsets for h, k, l directions ");
declareProperty(std::make_unique<PropertyWithValue<int>>("MaxOrder", 0, Direction::Input),
"Maximum order to apply ModVectors. Default = 0");
declareProperty("GetModVectorsFromUB", false, "If false Modulation Vectors will be read from input");
ModulationProperties::appendTo(this);
declareProperty(std::make_unique<PropertyWithValue<bool>>("CrossTerms", false, Direction::Input),
"Include cross terms (false)");
declareProperty("GetModVectorsFromUB", false, "If false Modulation Vectors will be read from input");
declareProperty("IncludeIntegerHKL", true, "If false order 0 peaks are not included in workspace (integer HKL)");
......@@ -81,7 +70,7 @@ void PredictSatellitePeaks::init() {
"Maximum wavelength limit at which to start looking for "
"single-crystal peaks.");
declareProperty(std::make_unique<PropertyWithValue<double>>("MinDSpacing", 0.1, Direction::Input),
"Minimum d-spacing of peaks to consider. Default = 1.0");
"Minimum d-spacing of peaks to consider. Default = 0.1");
declareProperty(std::make_unique<PropertyWithValue<double>>("MaxDSpacing", 100.0, Direction::Input),
"Maximum d-spacing of peaks to consider");
......@@ -109,11 +98,11 @@ void PredictSatellitePeaks::exec() {
return;
}
V3D offsets1 = getOffsetVector("ModVector1");
V3D offsets2 = getOffsetVector("ModVector2");
V3D offsets3 = getOffsetVector("ModVector3");
int maxOrder = getProperty("MaxOrder");
bool crossTerms = getProperty("CrossTerms");
V3D offsets1 = getOffsetVector(ModulationProperties::ModVector1);
V3D offsets2 = getOffsetVector(ModulationProperties::ModVector2);
V3D offsets3 = getOffsetVector(ModulationProperties::ModVector3);
int maxOrder = getProperty(ModulationProperties::MaxOrder);
bool crossTerms = getProperty(ModulationProperties::CrossTerms);
bool includeOrderZero = getProperty("IncludeIntegerHKL");
// boolean for only including order zero once
bool notOrderZero = false;
......@@ -140,7 +129,7 @@ void PredictSatellitePeaks::exec() {
const auto instrument = Peaks->getInstrument();
outPeaks = std::dynamic_pointer_cast<IPeaksWorkspace>(WorkspaceFactory::Instance().createPeaks(Peaks->id()));
outPeaks->setInstrument(instrument);
outPeaks->copyExperimentInfoFrom(Peaks.get());
outPeaks->mutableSample().setOrientedLattice(std::move(lattice));
Kernel::Matrix<double> goniometer;
......@@ -212,11 +201,11 @@ void PredictSatellitePeaks::exec() {
void PredictSatellitePeaks::exec_peaks() {
V3D offsets1 = getOffsetVector("ModVector1");
V3D offsets2 = getOffsetVector("ModVector2");
V3D offsets3 = getOffsetVector("ModVector3");
int maxOrder = getProperty("MaxOrder");
bool crossTerms = getProperty("CrossTerms");
V3D offsets1 = getOffsetVector(ModulationProperties::ModVector1);
V3D offsets2 = getOffsetVector(ModulationProperties::ModVector2);
V3D offsets3 = getOffsetVector(ModulationProperties::ModVector3);
int maxOrder = getProperty(ModulationProperties::MaxOrder);
bool crossTerms = getProperty(ModulationProperties::CrossTerms);
API::Sample sample = Peaks->mutableSample();
......@@ -251,7 +240,7 @@ void PredictSatellitePeaks::exec_peaks() {
const auto instrument = Peaks->getInstrument();
outPeaks = std::dynamic_pointer_cast<IPeaksWorkspace>(WorkspaceFactory::Instance().createPeaks(Peaks->id()));
outPeaks->setInstrument(instrument);
outPeaks->copyExperimentInfoFrom(Peaks.get());
outPeaks->mutableSample().setOrientedLattice(std::move(lattice));
vector<vector<int>> AlreadyDonePeaks;
......
......@@ -16,6 +16,7 @@
#include "MantidCrystal/CalculateUMatrix.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidKernel/V3D.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
......@@ -75,6 +76,47 @@ public:
AnalysisDataService::Instance().remove(WSName);
}
void test_exec_LeanElasticPeak() {
// Name of the output workspace.
std::string WSName("peaksCalculateUMatrix");
generatePeaks(WSName, true);
LeanElasticPeaksWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING(
ws = std::dynamic_pointer_cast<LeanElasticPeaksWorkspace>(AnalysisDataService::Instance().retrieve(WSName)));
TS_ASSERT(ws);
if (!ws)
return;
CalculateUMatrix alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("a", "2."));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("b", "3."));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("c", "4."));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("alpha", "90"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("beta", "90"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("gamma", "90"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("PeaksWorkspace", WSName));
TS_ASSERT_THROWS_NOTHING(alg.execute(););
TS_ASSERT(alg.isExecuted());
// Check that we set an oriented lattice
TS_ASSERT(ws->mutableSample().hasOrientedLattice());
// Check that the UB matrix is the same as in TOPAZ_3007.mat
OrientedLattice latt = ws->mutableSample().getOrientedLattice();
DblMatrix U(3, 3, false);
U[0][0] = sqrt(3.) * 0.5;
U[2][2] = sqrt(3.) * 0.5;
U[2][0] = 0.5;
U[0][2] = -0.5;
U[1][1] = 1.;
TS_ASSERT(latt.getU().equals(U, 1e-10));
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(WSName);
}
void test_fail() {
Mantid::Geometry::Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentRectangular2(1, 10);
inst->setName("SillyInstrument");
......@@ -167,7 +209,7 @@ private:
double ph(double H, double K, double L) { return atan2(-QYUB(H, K, L), -QXUB(H, K, L)); }
void generatePeaks(const std::string &WSName) {
void generatePeaks(const std::string &WSName, const bool leanElastic = false) {
setupUB();
double Hpeaks[9] = {0, 1, 1, 0, -1, -1, 1, -3, -2};
......@@ -183,8 +225,13 @@ private:
auto inst = ComponentCreationHelper::createCylInstrumentWithDetInGivenPositions(L2, theta, phi);
inst->setName("SillyInstrument");
auto pw = PeaksWorkspace_sptr(new PeaksWorkspace);
pw->setInstrument(inst);
IPeaksWorkspace_sptr pw;
if (leanElastic) {
pw = LeanElasticPeaksWorkspace_sptr(new LeanElasticPeaksWorkspace);
} else {
pw = PeaksWorkspace_sptr(new PeaksWorkspace);
pw->setInstrument(inst);
}
for (int i = 0; i <= 8; i++) {
Peak p(inst, i + 1, lambda[i], V3D(Hpeaks[i], Kpeaks[i], Lpeaks[i]));
pw->addPeak(p);
......
......@@ -7,6 +7,7 @@
#pragma once
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/IPeaksWorkspace.h"
#include "MantidAPI/Sample.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Timer.h"
......@@ -16,6 +17,8 @@
#include "MantidCrystal/LoadIsawUB.h"
#include "MantidCrystal/OptimizeLatticeForCellType.h"
#include "MantidDataHandling/LoadNexusProcessed.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
using namespace Mantid::Crystal;
......@@ -84,4 +87,63 @@ public:
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(WSName);
}
void test_exec_LeanElastic() {
// Name of the output workspace.
std::string WSName("peaks");
LoadNexusProcessed loader;
TS_ASSERT_THROWS_NOTHING(loader.initialize());
TS_ASSERT(loader.isInitialized());
loader.setPropertyValue("Filename", "TOPAZ_3007.peaks.nxs");
loader.setPropertyValue("OutputWorkspace", WSName);
TS_ASSERT(loader.execute());
TS_ASSERT(loader.isExecuted());
PeaksWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING(
ws = std::dynamic_pointer_cast<PeaksWorkspace>(AnalysisDataService::Instance().retrieve(WSName)));
TS_ASSERT(ws);
if (!ws)
return;
// Convert PeaksWorkspace to LeanElasticPeaksWorkspace
auto lpw = std::make_shared<LeanElasticPeaksWorkspace>();
for (auto peak : ws->getPeaks())
lpw->addPeak(peak);
AnalysisDataService::Instance().addOrReplace(WSName, lpw);
FindUBUsingFFT alg_fft;
TS_ASSERT_THROWS_NOTHING(alg_fft.initialize())
TS_ASSERT(alg_fft.isInitialized())
TS_ASSERT_THROWS_NOTHING(alg_fft.setPropertyValue("PeaksWorkspace", WSName));
TS_ASSERT_THROWS_NOTHING(alg_fft.setPropertyValue("MinD", "8.0"));
TS_ASSERT_THROWS_NOTHING(alg_fft.setPropertyValue("MaxD", "13.0"));
TS_ASSERT_THROWS_NOTHING(alg_fft.setPropertyValue("Tolerance", "0.15"));
TS_ASSERT_THROWS_NOTHING(alg_fft.execute(););
TS_ASSERT(alg_fft.isExecuted());
OptimizeLatticeForCellType alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("PeaksWorkspace", WSName));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("CellType", "Monoclinic"));
TS_ASSERT_THROWS_NOTHING(alg.execute(););
TS_ASSERT(alg.isExecuted());
// Check that we set an oriented lattice
TS_ASSERT(lpw->mutableSample().hasOrientedLattice());
// Check that the UB matrix is the same as in TOPAZ_3007.mat
OrientedLattice latt = lpw->mutableSample().getOrientedLattice();
double correct_UB[] = {-0.0451, 0.0406, -0.0125, 0.0016, -0.0031, 0.1158, 0.0574, 0.0322, 0.0275};
std::vector<double> UB_calculated = latt.getUB().getVector();
for (size_t i = 0; i < 9; i++) {
TS_ASSERT_DELTA(correct_UB[i], UB_calculated[i], 5e-4);
}
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(WSName);
}
};
......@@ -6,14 +6,14 @@
# SPDX - License - Identifier: GPL - 3.0 +
from mantid.api import (AlgorithmFactory, FileAction, FileProperty,
PythonAlgorithm, PropertyMode, ADSValidator,
WorkspaceGroup, WorkspaceProperty, MultipleExperimentInfos)
WorkspaceGroup, WorkspaceProperty, MultipleExperimentInfos,
IPeaksWorkspace)
from mantid.kernel import Direction, FloatBoundedValidator, StringListValidator, StringArrayProperty
from mantid.simpleapi import (DeleteWorkspace, IntegratePeaksMD,
SaveHKL, SaveReflections,
CreatePeaksWorkspace, CopySample,
AnalysisDataService,
CombinePeaksWorkspaces)
from mantid.dataobjects import PeaksWorkspace
AnalysisDataService, FilterPeaks,
CombinePeaksWorkspaces, mtd)
import numpy as np
......@@ -52,6 +52,9 @@ class HB3AIntegratePeaks(PythonAlgorithm):
self.declareProperty("ApplyLorentz", defaultValue=True,
doc="Whether the Lorentz correction should be applied to the integrated peaks")
self.declareProperty("RemoveZeroIntensity", defaultValue=True,
doc="If to remove peaks with 0 or less intensity from the output")
formats = StringListValidator()
formats.addAllowedValue("SHELX")
formats.addAllowedValue("Fullprof")
......@@ -92,8 +95,8 @@ class HB3AIntegratePeaks(PythonAlgorithm):
issues["InputWorkspace"] = "Workspace has the wrong number of dimensions"
for peak_ws in peak_workspaces:
if not isinstance(AnalysisDataService[peak_ws], PeaksWorkspace):
issues["PeaksWorkspace"] = "Workspace need to be a PeaksWorkspace"
if not isinstance(AnalysisDataService[peak_ws], IPeaksWorkspace):
issues["PeaksWorkspace"] = "Workspace need to be a PeaksWorkspace or LeanElasticPeaksWorkspace"
return issues
......@@ -105,6 +108,7 @@ class HB3AIntegratePeaks(PythonAlgorithm):
inner_radius = self.getProperty("BackgroundInnerRadius").value
outer_radius = self.getProperty("BackgroundOuterRadius").value
remove_0_intensity = self.getProperty("RemoveZeroIntensity").value
use_lorentz = self.getProperty("ApplyLorentz").value
multi_ws = len(input_workspaces) > 1
......@@ -129,7 +133,8 @@ class HB3AIntegratePeaks(PythonAlgorithm):
peaks_ws_name = output_workspace_name
CreatePeaksWorkspace(InstrumentWorkspace=input_workspaces[0],
NumberOfPeaks=0,
OutputWorkspace=peaks_ws_name)
OutputWorkspace=peaks_ws_name,
OutputType=mtd[peak_workspaces[0]].id().replace('sWorkspace',''))
CopySample(InputWorkspace=output_workspaces[0],
OutputWorkspace=peaks_ws_name,
CopyName=False,
......@@ -149,6 +154,10 @@ class HB3AIntegratePeaks(PythonAlgorithm):
lorentz = abs(np.sin(peak.getScattering() * np.cos(peak.getAzimuthal())))
peak.setIntensity(peak.getIntensity() * lorentz)
if remove_0_intensity:
FilterPeaks(InputWorkspace=peaks_ws_name, OutputWorkspace=peaks_ws_name,
FilterVariable='Intensity', FilterValue=0, Operator='>')
# Write output only if a file path was provided
if not self.getProperty("OutputFile").isDefault:
out_format = self.getProperty("OutputFormat").value
......
......@@ -5,8 +5,10 @@
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
from mantid.api import AlgorithmFactory, IMDWorkspaceProperty, IPeaksWorkspaceProperty, PythonAlgorithm, PropertyMode
from mantid.kernel import Direction, FloatPropertyWithValue, StringListValidator
from mantid.simpleapi import PredictPeaks, CloneMDWorkspace, CopySample, DeleteWorkspace, mtd
from mantid.kernel import Direction, FloatPropertyWithValue, StringListValidator, EnabledWhenProperty, PropertyCriterion
from mantid.simpleapi import (PredictPeaks, CloneMDWorkspace,
CopySample, DeleteWorkspace, PredictSatellitePeaks,
HFIRCalculateGoniometer, mtd)
import numpy as np
......@@ -16,7 +18,7 @@ class HB3APredictPeaks(PythonAlgorithm):
return "Crystal\\Peaks;Crystal\\UBMatrix"
def seeAlso(self):
return ["PredictPeaks", "HB3AFindPeaks", "HB3AIntegratePeaks", "HB3AAdjustSampleNorm"]
return ["PredictPeaks", "HB3AFindPeaks", "HB3AIntegratePeaks", "HB3AAdjustSampleNorm", "PredictSatellitePeaks"]
def name(self):
return "HB3APredictPeaks"
......@@ -51,6 +53,25 @@ class HB3APredictPeaks(PythonAlgorithm):
self.declareProperty("Wavelength", FloatPropertyWithValue.EMPTY_DBL,
doc="Wavelength value to use only if one was not found in the sample log")
self.declareProperty("SatellitePeaks", False, "If to Predict Satellite Peaks")
condition = EnabledWhenProperty("SatellitePeaks", PropertyCriterion.IsNotDefault)
self.declareProperty("ModVector1", "0.0,0.0,0.0", "Offsets for h, k, l directions")
self.declareProperty("ModVector2", "0.0,0.0,0.0", "Offsets for h, k, l directions")
self.declareProperty("ModVector3", "0.0,0.0,0.0", "Offsets for h, k, l directions")
self.declareProperty("MaxOrder", 0, "Maximum order to apply ModVectors. Default = 0")
self.declareProperty("GetModVectorsFromUB", False, "If false Modulation Vectors will be read from input")
self.declareProperty("CrossTerms", False, "Include cross terms (false)")
self.declareProperty("IncludeIntegerHKL", True, "If false order 0 peaks are not included in workspace (integer HKL)")
self.declareProperty("MinDSpacing", 1.0, "Minimum d-spacing of peaks to consider. Default = 1")
self.declareProperty("MaxDSpacing", 100.0, "Maximum d-spacing of peaks to consider. Default = 100")
self.setPropertySettings("ModVector1", condition)
self.setPropertySettings("ModVector2", condition)
self.setPropertySettings("ModVector3", condition)
self.setPropertySettings("MaxOrder", condition)
self.setPropertySettings("GetModVectorsFromUB", condition)
self.setPropertySettings("CrossTerms", condition)
self.setPropertySettings("IncludeIntegerHKL", condition)
self.declareProperty(IPeaksWorkspaceProperty("OutputWorkspace", defaultValue="", direction=Direction.Output,
optional=PropertyMode.Mandatory), doc="Output peaks workspace")
......@@ -87,41 +108,41 @@ class HB3APredictPeaks(PythonAlgorithm):
if input_ws.getNumExperimentInfo() == 0:
# Warn if we could extract a wavelength from the workspace
raise RuntimeWarning("No experiment info was found in input '{}'".format(input_ws.getName()))
else:
exp_info = input_ws.getExperimentInfo(0)
if exp_info.run().hasProperty("wavelength"):
wavelength = exp_info.run().getProperty("wavelength").value
if exp_info.run().hasProperty("omega") and exp_info.run().hasProperty("phi"):
gon = exp_info.run().getGoniometer().getEulerAngles('YZY')
if np.isclose(exp_info.run().getTimeAveragedStd("omega"), 0.0):
use_inner = True
min_angle = -exp_info.run().getLogData('phi').value.max()
max_angle = -exp_info.run().getLogData('phi').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
phi_log = -exp_info.run().getLogData('phi').value[0]
if np.isclose(phi_log + 180, gon[2]):
min_angle += 180
max_angle += 180
elif np.isclose(phi_log - 180, gon[2]):
min_angle -= 180
max_angle -= 180
elif np.isclose(exp_info.run().getTimeAveragedStd("phi"), 0.0):
use_inner = False
min_angle = -exp_info.run().getLogData('omega').value.max()
max_angle = -exp_info.run().getLogData('omega').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
omega_log = -exp_info.run().getLogData('omega').value[0]
if np.isclose(omega_log + 180, gon[0]):
min_angle += 180
max_angle += 180
elif np.isclose(omega_log - 180, gon[0]):
min_angle -= 180
max_angle -= 180
else:
self.log().warning("No appropriate goniometer rotation found, try anyway")
self.log().information("Using inner goniometer: {}".format(use_inner))
exp_info = input_ws.getExperimentInfo(0)
if exp_info.run().hasProperty("wavelength"):
wavelength = exp_info.run().getProperty("wavelength").value
if exp_info.run().hasProperty("omega") and exp_info.run().hasProperty("phi"):
gon = exp_info.run().getGoniometer().getEulerAngles('YZY')
if np.isclose(exp_info.run().getTimeAveragedStd("omega"), 0.0):
use_inner = True
min_angle = -exp_info.run().getLogData('phi').value.max()
max_angle = -exp_info.run().getLogData('phi').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
phi_log = -exp_info.run().getLogData('phi').value[0]
if np.isclose(phi_log + 180, gon[2]):
min_angle += 180
max_angle += 180
elif np.isclose(phi_log - 180, gon[2]):
min_angle -= 180
max_angle -= 180
elif np.isclose(exp_info.run().getTimeAveragedStd("phi"), 0.0):
use_inner = False
min_angle = -exp_info.run().getLogData('omega').value.max()
max_angle = -exp_info.run().getLogData('omega').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
omega_log = -exp_info.run().getLogData('omega').value[0]
if np.isclose(omega_log + 180, gon[0]):
min_angle += 180
max_angle += 180
elif np.isclose(omega_log - 180, gon[0]):
min_angle -= 180
max_angle -= 180
else:
self.log().warning("No appropriate goniometer rotation found, try anyway")
self.log().information("Using inner goniometer: {}".format(use_inner))
if not self.getProperty("Wavelength").isDefault:
wavelength = self.getProperty("Wavelength").value
......@@ -141,15 +162,38 @@ class HB3APredictPeaks(PythonAlgorithm):
CopyShape=False,
CopyLattice=True)
peaks = PredictPeaks(InputWorkspace=input_ws,
ReflectionCondition=reflection_condition,
CalculateGoniometerForCW=True,
Wavelength=wavelength,
FlipX=True,
InnerGoniometer=use_inner,
MinAngle=min_angle,
MaxAngle=max_angle,
OutputWorkspace=output_ws)
if self.getProperty("SatellitePeaks").value: