Commit 8de22e32 authored by Matthew Bowles's avatar Matthew Bowles
Browse files

initial fit for single/multi ion/spectra Re #20358

parent c3d1dd66
......@@ -1711,7 +1711,7 @@ void Algorithm::execMasterOnly() {
* non-master ranks in master-only execution. */
void Algorithm::execNonMaster() {
// If there is no output we can simply do nothing.
if (m_pureOutputWorkspaceProps.size() == 0)
if (m_pureOutputWorkspaceProps.empty())
return;
// Does Algorithm have exactly one input and one output workspace property?
if (m_inputWorkspaceProps.size() == 1 &&
......
......@@ -148,7 +148,7 @@ DetectorSearcher::searchUsingNearestNeighbours(const V3D &q) {
// find where this Q vector should intersect with "extended" space
const auto neighbours =
m_detectorCacheSearch->findNearest(Eigen::Vector3d(q[0], q[1], q[2]), 5);
if (neighbours.size() == 0)
if (neighbours.empty())
return std::make_tuple(false, 0);
const auto result = checkInteceptWithNeighbours(detectorDir, neighbours);
......
......@@ -10,7 +10,7 @@ IndexTypeProperty::IndexTypeProperty(const std::string &name,
if (indexType & IndexType::SpectrumNum)
m_allowedValues.push_back("SpectrumNumber");
if (m_allowedValues.size() == 0)
if (m_allowedValues.empty())
throw std::invalid_argument("Argument indexType incorrectly specified");
m_value = m_allowedValues[0];
......
......@@ -239,7 +239,7 @@ void MatrixWorkspace::initialize(const std::size_t &NVectors,
void MatrixWorkspace::initialize(const std::size_t &NVectors,
const HistogramData::Histogram &histogram) {
// Check validity of arguments
if (NVectors == 0 || histogram.x().size() == 0) {
if (NVectors == 0 || histogram.x().empty()) {
throw std::out_of_range(
"All arguments to init must be positive and non-zero");
}
......
......@@ -8,6 +8,9 @@ namespace Mantid {
namespace API {
class SpectrumInfo;
}
namespace Geometry {
class ReferenceFrame;
}
namespace HistogramData {
class HistogramX;
class HistogramY;
......@@ -88,6 +91,8 @@ private:
// convert to momentum transfer
Mantid::API::MatrixWorkspace_sptr
convertToQ(Mantid::API::MatrixWorkspace_sptr inputWS);
// Get the twoTheta width of a given detector
double getDetectorTwoThetaRange(const size_t spectrumIdx);
// Utility function to create name for diagnostic workspaces
std::string createDebugWorkspaceName(const std::string &inputName);
// Utility function to output a diagnostic workspace to the ADS
......@@ -152,6 +157,7 @@ private:
API::MatrixWorkspace_sptr m_runWS;
const API::SpectrumInfo *m_spectrumInfo;
boost::shared_ptr<const Mantid::Geometry::ReferenceFrame> m_refFrame;
bool m_convertUnits; // convert the input workspace to lambda
bool m_normaliseMonitors; // normalise by monitors and direct beam
bool m_normaliseTransmission; // transmission or algorithmic correction
......
......@@ -228,7 +228,7 @@ void AddSampleLog::addTimeSeriesProperty(Run &run_obj,
is_int_series = true;
} else if (prop_number_type == autoTypeOption) {
// auto type. by default
if (prop_value.size() == 0)
if (prop_value.empty())
g_log.warning("For sample log in TimeSeriesProperty and values are given "
"by MarixWorkspace, the default data type "
"is double.");
......@@ -248,8 +248,8 @@ void AddSampleLog::addTimeSeriesProperty(Run &run_obj,
// check using workspace or some specified start value
std::string tsp_ws_name = getPropertyValue("TimeSeriesWorkspace");
bool use_ws = tsp_ws_name.size() > 0;
bool use_single_value = prop_value.size() > 0;
bool use_ws = !tsp_ws_name.empty();
bool use_single_value = !prop_value.empty();
if (use_ws && use_single_value) {
throw std::runtime_error("Both TimeSeries workspace and sing value are "
"specified. It is not allowed.");
......
......@@ -331,7 +331,7 @@ bool ChangeTimeZero::checkForDateTime(const std::string &val) const {
// Hedge for bad lexical casts in the DateTimeValidator
try {
DateTimeValidator validator = DateTimeValidator();
isDateTime = validator.isValid(val) == "";
isDateTime = validator.isValid(val).empty();
} catch (...) {
isDateTime = false;
}
......
......@@ -105,7 +105,7 @@ void ConvertAxisByFormula::exec() {
RefAxis *refAxisPtr = dynamic_cast<RefAxis *>(axisPtr);
if (refAxisPtr != nullptr) {
CommonBinsValidator sameBins;
if (sameBins.isValid(outputWs) != "") {
if (!sameBins.isValid(outputWs).empty()) {
isRaggedBins = true;
}
isRefAxis = true;
......@@ -245,14 +245,14 @@ void ConvertAxisByFormula::exec() {
}
// Set the Unit of the Axis
if ((axisUnits != "") || (axisTitle != "")) {
if ((!axisUnits.empty()) || (!axisTitle.empty())) {
try {
axisPtr->unit() = UnitFactory::Instance().create(axisUnits);
} catch (Exception::NotFoundError &) {
if (axisTitle == "") {
if (axisTitle.empty()) {
axisTitle = axisPtr->unit()->caption();
}
if (axisUnits == "") {
if (axisUnits.empty()) {
axisUnits = axisPtr->unit()->label();
}
axisPtr->unit() = boost::make_shared<Units::Label>(axisTitle, axisUnits);
......
......@@ -337,7 +337,7 @@ ConvertUnits::convertQuickly(API::MatrixWorkspace_const_sptr inputWS,
// First a quick check using the validator
CommonBinsValidator sameBins;
bool commonBoundaries = false;
if (sameBins.isValid(inputWS) == "") {
if (sameBins.isValid(inputWS).empty()) {
commonBoundaries = WorkspaceHelpers::commonBoundaries(*inputWS);
// Only do the full check if the quick one passes
if (commonBoundaries) {
......
......@@ -299,7 +299,7 @@ void DiffractionEventCalibrateDetectors::exec() {
std::vector<boost::shared_ptr<RectangularDetector>> detList;
// --------- Loading only one bank ----------------------------------
std::string onebank = getProperty("BankName");
bool doOneBank = (onebank != "");
bool doOneBank = (!onebank.empty());
for (int i = 0; i < inst->nelements(); i++) {
boost::shared_ptr<RectangularDetector> det;
boost::shared_ptr<ICompAssembly> assem;
......
......@@ -109,7 +109,7 @@ void DiffractionFocussing2::exec() {
throw std::invalid_argument("Workspace Invalid Spacing/UnitID");
}
// --- Do we need to read the grouping workspace? ----
if (groupingFileName != "") {
if (!groupingFileName.empty()) {
progress(0.01, "Reading grouping file");
IAlgorithm_sptr childAlg = createChildAlgorithm("CreateGroupingWorkspace");
childAlg->setProperty(
......
......@@ -429,7 +429,7 @@ void ExportTimeSeriesLog::calculateFirstDerivative(bool is_event_ws) {
// error message
std::string errmsg = errmsg_ss.str();
if (errmsg.size() > 0)
if (!errmsg.empty())
g_log.error(errmsg);
return;
......
......@@ -77,12 +77,12 @@ void FilterByTime::exec() {
start_str = getPropertyValue("AbsoluteStartTime");
stop_str = getPropertyValue("AbsoluteStopTime");
if ((start_str != "") && (stop_str != "") && (start_dbl <= 0.0) &&
if ((!start_str.empty()) && (!stop_str.empty()) && (start_dbl <= 0.0) &&
(stop_dbl <= 0.0)) {
// Use the absolute string
start = DateAndTime(start_str);
stop = DateAndTime(stop_str);
} else if ((start_str == "") && (stop_str == "") &&
} else if ((start_str.empty()) && (stop_str.empty()) &&
((start_dbl > 0.0) || (stop_dbl > 0.0))) {
// Use the relative times in seconds.
DateAndTime first = inputWS->getFirstPulseTime();
......
......@@ -73,12 +73,12 @@ void FilterByTime2::exec() {
std::string absstoptime = this->getProperty("AbsoluteStopTime");
std::string start, stop;
if ((absstarttime != "") && (absstoptime != "") && (starttime <= 0.0) &&
if ((!absstarttime.empty()) && (!absstoptime.empty()) && (starttime <= 0.0) &&
(stoptime <= 0.0)) {
// Use the absolute string
start = absstarttime;
stop = absstoptime;
} else if ((absstarttime != "" || absstoptime != "") &&
} else if ((!absstarttime.empty() || !absstoptime.empty()) &&
(starttime > 0.0 || stoptime > 0.0)) {
throw std::invalid_argument(
"It is not allowed to provide both absolute time and relative time.");
......
......@@ -792,7 +792,7 @@ void FilterEvents::convertSplittersWorkspaceToVectors() {
Kernel::SplittingInterval splitter = m_splitters[i_splitter];
int64_t start_time_i64 = splitter.start().totalNanoseconds();
int64_t stop_time_i64 = splitter.stop().totalNanoseconds();
if (m_vecSplitterTime.size() == 0) {
if (m_vecSplitterTime.empty()) {
// first entry: add
m_vecSplitterTime.push_back(start_time_i64);
m_vecSplitterTime.push_back(stop_time_i64);
......@@ -949,7 +949,7 @@ void FilterEvents::processTableSplittersWorkspace() {
int64_t stop_64 =
filter_shift_time + static_cast<int64_t>(stop_time * 1.E9);
if (m_vecSplitterTime.size() == 0) {
if (m_vecSplitterTime.empty()) {
// first splitter: push the start time to vector
m_vecSplitterTime.push_back(start_64);
} else if (start_64 - m_vecSplitterTime.back() > TOLERANCE) {
......@@ -973,7 +973,7 @@ void FilterEvents::processTableSplittersWorkspace() {
// convert string-target to integer target
bool addnew = false;
int int_target(-1);
if (m_targetIndexMap.size() == 0) {
if (m_targetIndexMap.empty()) {
addnew = true;
} else {
std::map<std::string, int>::iterator mapiter =
......
......@@ -80,8 +80,8 @@ void GenerateIPythonNotebook::exec() {
}
// Need at least a start time to do time filter
if (startTime != "") {
if (endTime == "") {
if (!startTime.empty()) {
if (endTime.empty()) {
// If no end time was given then filter up to now
view->filterBetweenExecDate(DateAndTime(startTime));
} else {
......
......@@ -79,8 +79,8 @@ void GeneratePythonScript::exec() {
}
// Need at least a start time to do time filter
if (startTime != "") {
if (endTime == "") {
if (!startTime.empty()) {
if (endTime.empty()) {
// If no end time was given then filter up to now
view->filterBetweenExecDate(DateAndTime(startTime));
} else {
......
......@@ -213,7 +213,7 @@ void Q1DWeighted::exec() {
// For reference - in the case where we don't use sub-pixels, simply
// use:
// double sinTheta = sin( spectrumInfo.twoTheta(i)/2.0 );
V3D pos = spectrumInfo.position(i) - V3D(sub_x, sub_y, 0.0);
V3D pos = spectrumInfo.position(i) - V3D(sub_x, sub_y, 0.0) - samplePos;
double sinTheta = sin(0.5 * pos.angle(beamLine));
double factor = fmp * sinTheta;
double q = factor * 2.0 / (XIn[j] + XIn[j + 1]);
......
......@@ -4,18 +4,22 @@
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/Objects/BoundingBox.h"
#include "MantidHistogramData/LinearGenerator.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/Unit.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include <algorithm>
#include <boost/lexical_cast.hpp>
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::HistogramData;
using namespace Mantid::Indexing;
......@@ -37,27 +41,6 @@ double getDetectorTwoTheta(const SpectrumInfo *spectrumInfo,
return spectrumInfo->signedTwoTheta(spectrumIdx);
}
/** Get the twoTheta angle range for the top/bottom of the detector associated
* with the given spectrum
*
* @param spectrumInfo : the spectrum info
* @param spectrumIdx : the workspace index of the spectrum
* @return : the twoTheta angle in radians
*/
double getDetectorTwoThetaRange(const SpectrumInfo *spectrumInfo,
const size_t spectrumIdx) {
// Assume the range covered by this pixel is the diff between this
// pixel's twoTheta and the next/prev pixel)
double twoTheta = getDetectorTwoTheta(spectrumInfo, spectrumIdx);
double bTwoTheta = 0;
if (spectrumIdx + 1 < spectrumInfo->size()) {
bTwoTheta = getDetectorTwoTheta(spectrumInfo, spectrumIdx + 1) - twoTheta;
}
return bTwoTheta;
}
/** Get the start/end of the lambda range for the detector associated
* with the given spectrum
*
......@@ -328,6 +311,34 @@ std::string createProcessingCommandsFromDetectorWS(
return result;
}
/**
* Get the topbottom extent of a detector for the given axis
*
* @param axis [in] : the axis to get the extent for
* @param top [in] : if true, get the max extent, or min otherwise
* @return : the max/min extent on the given axis
*/
double getBoundingBoxExtent(const BoundingBox &boundingBox,
const PointingAlong axis, const bool top) {
double result = 0.0;
switch (axis) {
case X:
result = top ? boundingBox.xMax() : boundingBox.xMin();
break;
case Y:
result = top ? boundingBox.yMax() : boundingBox.yMin();
break;
case Z:
result = top ? boundingBox.zMax() : boundingBox.zMin();
break;
default:
throw std::runtime_error("Axis is not X/Y/Z");
break;
}
return result;
}
}
// Register the algorithm into the AlgorithmFactory
......@@ -434,7 +445,10 @@ void ReflectometryReductionOne2::exec() {
throw std::invalid_argument(
"InputWorkspace must have units of TOF or Wavelength");
// Cache the spectrum info and reference frame
m_spectrumInfo = &m_runWS->spectrumInfo();
auto instrument = m_runWS->getInstrument();
m_refFrame = instrument->getReferenceFrame();
// Find and cache detector groups and theta0
findDetectorGroups();
......@@ -463,6 +477,44 @@ void ReflectometryReductionOne2::exec() {
setProperty("OutputWorkspace", IvsQ);
}
/** Get the twoTheta angle range for the top/bottom of the detector associated
* with the given spectrum
*
* @param spectrumIdx : the workspace index of the spectrum
* @return : the twoTheta range in radians
*/
double
ReflectometryReductionOne2::getDetectorTwoThetaRange(const size_t spectrumIdx) {
double bTwoTheta = 0;
// Get the sample->detector distance along the beam
const V3D detSample =
m_spectrumInfo->position(spectrumIdx) - m_spectrumInfo->samplePosition();
const double beamOffset =
m_refFrame->vecPointingAlongBeam().scalar_prod(detSample);
// Get the bounding box for this detector/group
BoundingBox boundingBox;
auto detector = m_runWS->getDetector(spectrumIdx);
detector->getBoundingBox(boundingBox);
// Get the top and bottom on the axis pointing up
const double top =
getBoundingBoxExtent(boundingBox, m_refFrame->pointingUp(), true);
const double bottom =
getBoundingBoxExtent(boundingBox, m_refFrame->pointingUp(), false);
// Calculate the difference in twoTheta between the top and bottom
const double twoThetaTop = std::atan(top / beamOffset);
const double twoThetaBottom = std::atan(bottom / beamOffset);
bTwoTheta = twoThetaTop - twoThetaBottom;
// We must have non-zero width to project a range
if (bTwoTheta < Tolerance) {
throw std::runtime_error("Error calculating pixel size.");
}
return bTwoTheta;
}
/**
* Utility function to create a unique workspace name for diagnostic outputs
* based on a given input workspace name
......@@ -856,9 +908,19 @@ void ReflectometryReductionOne2::findDetectorGroups() {
return a.front() < b.front();
});
if (m_detectorGroups.size() == 0) {
if (m_detectorGroups.empty()) {
throw std::runtime_error("Invalid processing instructions");
}
for (const auto &group : m_detectorGroups) {
for (const auto &spIdx : group) {
if (spIdx > m_spectrumInfo->size() - 1) {
throw std::runtime_error(
"ProcessingInstructions contains an out-of-range index: " +
std::to_string(spIdx));
}
}
}
}
/**
......@@ -993,8 +1055,7 @@ void ReflectometryReductionOne2::findIvsLamRange(
const size_t spIdxMin = detectors.front();
const double twoThetaMin = getDetectorTwoTheta(m_spectrumInfo, spIdxMin);
const double bTwoThetaMin =
getDetectorTwoThetaRange(m_spectrumInfo, spIdxMin);
const double bTwoThetaMin = getDetectorTwoThetaRange(spIdxMin);
// For bLambda, use the average bin size for this spectrum
auto xValues = detectorWS->x(spIdxMin);
double bLambda = (xValues[xValues.size() - 1] - xValues[0]) /
......@@ -1004,8 +1065,7 @@ void ReflectometryReductionOne2::findIvsLamRange(
const size_t spIdxMax = detectors.back();
const double twoThetaMax = getDetectorTwoTheta(m_spectrumInfo, spIdxMax);
const double bTwoThetaMax =
getDetectorTwoThetaRange(m_spectrumInfo, spIdxMax);
const double bTwoThetaMax = getDetectorTwoThetaRange(spIdxMax);
xValues = detectorWS->x(spIdxMax);
bLambda = (xValues[xValues.size() - 1] - xValues[0]) /
static_cast<int>(xValues.size());
......@@ -1092,7 +1152,7 @@ ReflectometryReductionOne2::sumInQ(MatrixWorkspace_sptr detectorWS) {
for (auto spIdx : detectors) {
// Get the angle of this detector and its size in twoTheta
const double twoTheta = getDetectorTwoTheta(m_spectrumInfo, spIdx);
const double bTwoTheta = getDetectorTwoThetaRange(m_spectrumInfo, spIdx);
const double bTwoTheta = getDetectorTwoThetaRange(spIdx);
// Check X length is Y length + 1
const auto &inputX = detectorWS->x(spIdx);
......@@ -1228,11 +1288,19 @@ void ReflectometryReductionOne2::sumInQShareCounts(
}
// Add a share of the input counts to this bin based on the proportion of
// overlap.
if (totalWidth > Tolerance) {
// Share counts out proportionally based on the overlap of this range
const double overlapWidth =
std::min({bLambda, lambdaMax - binStart, binEnd - lambdaMin});
const double fraction = overlapWidth / totalWidth;
outputY[outIdx] += inputCounts * fraction;
outputE[outIdx] += inputErr * fraction;
} else {
// Projection to a single value. Put all counts in the overlapping output
// bin.
outputY[outIdx] += inputCounts;
outputE[outIdx] += inputCounts;
}
}
}
......
......@@ -61,15 +61,15 @@ std::map<std::string, std::string> RenameWorkspaces::validateInputs() {
std::string suffix = getPropertyValue("Suffix");
// Check properties
if (newWsName.empty() && prefix == "" && suffix == "") {
if (newWsName.empty() && prefix.empty() && suffix.empty()) {
errorList["WorkspaceNames"] =
"No list of Workspace names, prefix or suffix has been supplied.";
}
if (!newWsName.empty() && (prefix != "" || suffix != "")) {
if (!newWsName.empty() && (!prefix.empty() || !suffix.empty())) {
errorList["WorkspaceNames"] = "Both a list of workspace names and a prefix "
"or suffix has been supplied.";
if (prefix != "") {
if (!prefix.empty()) {
errorList["Prefix"] = "Both a list of workspace names and a prefix "
"or suffix has been supplied.";
} else {
......
Supports Markdown
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