#include "MantidAlgorithms/ReflectometryReductionOneAuto2.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAlgorithms/BoostOptionalToAlgorithmProperty.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/make_unique.h"
namespace Mantid {
namespace Algorithms {
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ReflectometryReductionOneAuto2)
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string ReflectometryReductionOneAuto2::name() const {
return "ReflectometryReductionOneAuto";
}
/// Algorithm's version for identification. @see Algorithm::version
int ReflectometryReductionOneAuto2::version() const { return 2; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string ReflectometryReductionOneAuto2::category() const {
return "Reflectometry\\ISIS";
}
/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string ReflectometryReductionOneAuto2::summary() const {
return "Reduces a single TOF/Lambda reflectometry run into a mod Q vs I/I0 "
"workspace attempting to pick instrument parameters for missing "
"properties";
}
/** Validate transmission runs
*
* @return :: result of the validation as a map
*/
std::map<std::string, std::string>
ReflectometryReductionOneAuto2::validateInputs() {
std::map<std::string, std::string> results;
// Validate transmission runs only if our input workspace is a group
if (!checkGroups())
return results;
auto group = AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(
getPropertyValue("InputWorkspace"));
if (!group)
return results;
// First transmission run
const std::string firstStr = getPropertyValue("FirstTransmissionRun");
if (!firstStr.empty()) {
auto firstTransGroup =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(firstStr);
// If it is not a group, we don't need to validate its size
if (!firstTransGroup)
return results;
const bool polarizationCorrections =
getPropertyValue("PolarizationAnalysis") != "None";
if (group->size() != firstTransGroup->size() && !polarizationCorrections) {
// If they are not the same size then we cannot associate a transmission
// group member with every input group member.
results["FirstTransmissionRun"] =
"FirstTransmissionRun group must be the "
"same size as the InputWorkspace group "
"when polarization analysis is 'None'.";
}
}
// The same for the second transmission run
const std::string secondStr = getPropertyValue("SecondTransmissionRun");
if (!secondStr.empty()) {
auto secondTransGroup =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(secondStr);
// If it is not a group, we don't need to validate its size
if (!secondTransGroup)
return results;
const bool polarizationCorrections =
getPropertyValue("PolarizationAnalysis") != "None";
if (group->size() != secondTransGroup->size() && !polarizationCorrections) {
results["SecondTransmissionRun"] =
"SecondTransmissionRun group must be the "
"same size as the InputWorkspace group "
"when polarization analysis is 'None'.";
}
}
return results;
}
/** Initialize the algorithm's properties.
*/
void ReflectometryReductionOneAuto2::init() {
// Input ws
declareProperty(
make_unique<WorkspaceProperty<MatrixWorkspace>>(
"InputWorkspace", "", Direction::Input, PropertyMode::Mandatory),
"Input run in TOF or wavelength");
// Reduction type
initReductionProperties();
// Analysis mode
const std::vector<std::string> analysisMode{"PointDetectorAnalysis",
"MultiDetectorAnalysis"};
auto analysisModeValidator =
boost::make_shared<StringListValidator>(analysisMode);
declareProperty("AnalysisMode", analysisMode[0], analysisModeValidator,
"Analysis mode. This property is only used when "
"ProcessingInstructions is not set.",
Direction::Input);
// Processing instructions
declareProperty(make_unique<PropertyWithValue<std::string>>(
"ProcessingInstructions", "", Direction::Input),
"Grouping pattern of workspace indices to yield only the"
" detectors of interest. See GroupDetectors for syntax.");
// Theta
declareProperty("ThetaIn", Mantid::EMPTY_DBL(), "Angle in degrees",
Direction::Input);
// ThetaLogName
declareProperty("ThetaLogName", "",
"The name ThetaIn can be found in the run log as");
// Whether to correct detectors
declareProperty(
make_unique<PropertyWithValue<bool>>("CorrectDetectors", true,
Direction::Input),
"Moves detectors to twoTheta if ThetaIn or ThetaLogName is given");
// Detector position correction type
const std::vector<std::string> correctionType{"VerticalShift",
"RotateAroundSample"};
auto correctionTypeValidator = boost::make_shared<CompositeValidator>();
correctionTypeValidator->add(
boost::make_shared<MandatoryValidator<std::string>>());
correctionTypeValidator->add(
boost::make_shared<StringListValidator>(correctionType));
declareProperty(
"DetectorCorrectionType", correctionType[0], correctionTypeValidator,
"When correcting detector positions, this determines whether detectors"
"should be shifted vertically or rotated around the sample position.",
Direction::Input);
setPropertySettings("DetectorCorrectionType",
make_unique<Kernel::EnabledWhenProperty>(
"CorrectDetectors", IS_EQUAL_TO, "1"));
// Wavelength limits
declareProperty("WavelengthMin", Mantid::EMPTY_DBL(),
"Wavelength Min in angstroms", Direction::Input);
declareProperty("WavelengthMax", Mantid::EMPTY_DBL(),
"Wavelength Max in angstroms", Direction::Input);
// Monitor properties
initMonitorProperties();
// Normalization by integrated monitors
declareProperty("NormalizeByIntegratedMonitors", true,
"Normalize by dividing by the integrated monitors.");
// Init properties for transmission normalization
initTransmissionProperties();
// Init properties for algorithmic corrections
initAlgorithmicProperties(true);
// Momentum transfer properties
initMomentumTransferProperties();
// Polarization correction
std::vector<std::string> propOptions = {"None", "PA", "PNR", "ParameterFile"};
declareProperty("PolarizationAnalysis", "None",
boost::make_shared<StringListValidator>(propOptions),
"Polarization analysis mode.");
declareProperty(
Kernel::make_unique<ArrayProperty<double>>("Pp", Direction::Input),
"Effective polarizing power of the polarizing system. "
"Expressed as a ratio 0 < Pp < 1");
declareProperty(
Kernel::make_unique<ArrayProperty<double>>("Ap", Direction::Input),
"Effective polarizing power of the analyzing system. "
"Expressed as a ratio 0 < Ap < 1");
declareProperty(
Kernel::make_unique<ArrayProperty<double>>("Rho", Direction::Input),
"Ratio of efficiencies of polarizer spin-down to polarizer "
"spin-up. This is characteristic of the polarizer flipper. "
"Values are constants for each term in a polynomial "
"expression.");
declareProperty(
Kernel::make_unique<ArrayProperty<double>>("Alpha", Direction::Input),
"Ratio of efficiencies of analyzer spin-down to analyzer "
"spin-up. This is characteristic of the analyzer flipper. "
"Values are factors for each term in a polynomial "
"expression.");
setPropertyGroup("PolarizationAnalysis", "Polarization Corrections");
setPropertyGroup("Pp", "Polarization Corrections");
setPropertyGroup("Ap", "Polarization Corrections");
setPropertyGroup("Rho", "Polarization Corrections");
setPropertyGroup("Alpha", "Polarization Corrections");
// Init properties for diagnostics
initDebugProperties();
// Output workspace in Q
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspaceBinned", "", Direction::Output),
"Output workspace in Q (rebinned workspace)");
// Output workspace in Q (unbinned)
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"Output workspace in Q (native binning)");
// Output workspace in wavelength
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspaceWavelength", "", Direction::Output),
"Output workspace in wavelength");
}
/** Execute the algorithm.
*/
void ReflectometryReductionOneAuto2::exec() {
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
auto instrument = inputWS->getInstrument();
IAlgorithm_sptr alg = createChildAlgorithm("ReflectometryReductionOne");
alg->initialize();
// Mandatory properties
alg->setProperty("SummationType", getPropertyValue("SummationType"));
alg->setProperty("ReductionType", getPropertyValue("ReductionType"));
alg->setProperty("IncludePartialBins",
getPropertyValue("IncludePartialBins"));
alg->setProperty("Diagnostics", getPropertyValue("Diagnostics"));
double wavMin = checkForMandatoryInstrumentDefault<double>(
this, "WavelengthMin", instrument, "LambdaMin");
alg->setProperty("WavelengthMin", wavMin);
double wavMax = checkForMandatoryInstrumentDefault<double>(
this, "WavelengthMax", instrument, "LambdaMax");
alg->setProperty("WavelengthMax", wavMax);
const auto instructions =
populateProcessingInstructions(alg, instrument, inputWS);
// Now that we know the detectors of interest, we can move them if necessary
// (i.e. if theta is given). If not, we calculate theta from the current
// detector positions
bool correctDetectors = getProperty("CorrectDetectors");
double theta;
if (!getPointerToProperty("ThetaIn")->isDefault()) {
theta = getProperty("ThetaIn");
} else if (!getPropertyValue("ThetaLogName").empty()) {
theta = getThetaFromLogs(inputWS, getPropertyValue("ThetaLogName"));
} else {
// Calculate theta from detector positions
theta = calculateTheta(instructions, inputWS);
// Never correct detector positions if ThetaIn or ThetaLogName is not
// specified
correctDetectors = false;
}
// Pass theta to the child algorithm
alg->setProperty("ThetaIn", theta);
if (correctDetectors) {
inputWS = correctDetectorPositions(instructions, inputWS, 2 * theta);
}
// Optional properties
populateMonitorProperties(alg, instrument);
alg->setPropertyValue("NormalizeByIntegratedMonitors",
getPropertyValue("NormalizeByIntegratedMonitors"));
bool transRunsFound = populateTransmissionProperties(alg);
if (transRunsFound)
alg->setProperty("StrictSpectrumChecking",
getPropertyValue("StrictSpectrumChecking"));
else
populateAlgorithmicCorrectionProperties(alg, instrument);
alg->setProperty("InputWorkspace", inputWS);
alg->execute();
MatrixWorkspace_sptr IvsLam = alg->getProperty("OutputWorkspaceWavelength");
MatrixWorkspace_sptr IvsQ = alg->getProperty("OutputWorkspace");
std::vector<double> params;
MatrixWorkspace_sptr IvsQB = rebinAndScale(IvsQ, theta, params);
setProperty("OutputWorkspaceWavelength", IvsLam);
setProperty("OutputWorkspace", IvsQ);
setProperty("OutputWorkspaceBinned", IvsQB);
// Set other properties so they can be updated in the Reflectometry interface
setProperty("ThetaIn", theta);
if (!params.empty()) {
if (params.size() == 3) {
setProperty("MomentumTransferMin", params[0]);
setProperty("MomentumTransferStep", -params[1]);
setProperty("MomentumTransferMax", params[2]);
} else {
setProperty("MomentumTransferMin", IvsQ->x(0).front());
setProperty("MomentumTransferMax", IvsQ->x(0).back());
setProperty("MomentumTransferStep", -params[0]);
}
}
if (getPointerToProperty("ScaleFactor")->isDefault())
setProperty("ScaleFactor", 1.0);
}
/** Returns the detectors of interest, specified via processing instructions.
* Note that this returns the names of the parent detectors of the first and
* last spectrum indices in the processing instructions. It is assumed that all
* the interim detectors have the same parent.
*
* @param instructions :: processing instructions defining detectors of interest
* @param inputWS :: the input workspace
* @return :: the names of the detectors of interest
*/
std::vector<std::string> ReflectometryReductionOneAuto2::getDetectorNames(
const std::string &instructions, MatrixWorkspace_sptr inputWS) {
std::vector<std::string> wsIndices;
boost::split(wsIndices, instructions, boost::is_any_of(":,-+"));
// vector of comopnents
std::vector<std::string> detectors;
try {
for (const auto &wsIndex : wsIndices) {
size_t index = boost::lexical_cast<size_t>(wsIndex);
auto detector = inputWS->getDetector(index);
auto parent = detector->getParent();
if (parent) {
auto parentType = parent->type();
auto detectorName = (parentType == "Instrument") ? detector->getName()
: parent->getName();
detectors.push_back(detectorName);
}
}
} catch (boost::bad_lexical_cast &) {
throw std::runtime_error("Invalid processing instructions: " +
instructions);
}
return detectors;
}
/** Correct an instrument component by shifting it vertically or
* rotating it around the sample.
*
* @param instructions :: processing instructions defining the detectors of
* interest
* @param inputWS :: the input workspace
* @param twoTheta :: the angle to move detectors to
* @return :: the corrected workspace
*/
MatrixWorkspace_sptr ReflectometryReductionOneAuto2::correctDetectorPositions(
const std::string &instructions, MatrixWorkspace_sptr inputWS,
const double twoTheta) {
auto detectorsOfInterest = getDetectorNames(instructions, inputWS);
// Detectors of interest may be empty. This happens for instance when we input
// a workspace that was previously reduced using this algorithm. In this case,
// we shouldn't correct the detector positions
if (detectorsOfInterest.empty())
return inputWS;
const std::set<std::string> detectorSet(detectorsOfInterest.begin(),
detectorsOfInterest.end());
const std::string correctionType = getProperty("DetectorCorrectionType");
MatrixWorkspace_sptr corrected = inputWS;
for (const auto &detector : detectorSet) {
IAlgorithm_sptr alg =
createChildAlgorithm("SpecularReflectionPositionCorrect");
alg->setProperty("InputWorkspace", corrected);
alg->setProperty("TwoTheta", twoTheta);
alg->setProperty("DetectorCorrectionType", correctionType);
alg->setProperty("DetectorComponentName", detector);
alg->execute();
corrected = alg->getProperty("OutputWorkspace");
}
return corrected;
}
/** Calculate the theta value of the detector of interest specified via
* processing instructions
*
* @param instructions :: processing instructions defining the detectors of
* interest
* @param inputWS :: the input workspace
* @return :: the angle of the detector (only the first detector is considered)
*/
double
ReflectometryReductionOneAuto2::calculateTheta(const std::string &instructions,
MatrixWorkspace_sptr inputWS) {
const auto detectorsOfInterest = getDetectorNames(instructions, inputWS);
// Detectors of interest may be empty. This happens for instance when we input
// a workspace that was previously reduced using this algorithm. In this case,
// we can't calculate theta
if (detectorsOfInterest.empty())
return 0.0;
IAlgorithm_sptr alg =
createChildAlgorithm("SpecularReflectionCalculateTheta");
alg->setProperty("InputWorkspace", inputWS);
alg->setProperty("DetectorComponentName", detectorsOfInterest[0]);
alg->execute();
const double theta = alg->getProperty("TwoTheta");
// Take a factor of 0.5 of the detector position, which is expected to be at
// 2 * theta
return theta * 0.5;
}
/** Set algorithmic correction properties
*
* @param alg :: ReflectometryReductionOne algorithm
* @param instrument :: The instrument attached to the workspace
*/
void ReflectometryReductionOneAuto2::populateAlgorithmicCorrectionProperties(
IAlgorithm_sptr alg, Instrument_const_sptr instrument) {
// With algorithmic corrections, monitors should not be integrated, see below
const std::string correctionAlgorithm = getProperty("CorrectionAlgorithm");
if (correctionAlgorithm == "PolynomialCorrection") {
alg->setProperty("NormalizeByIntegratedMonitors", false);
alg->setProperty("CorrectionAlgorithm", "PolynomialCorrection");
alg->setPropertyValue("Polynomial", getPropertyValue("Polynomial"));
} else if (correctionAlgorithm == "ExponentialCorrection") {
alg->setProperty("NormalizeByIntegratedMonitors", false);
alg->setProperty("CorrectionAlgorithm", "ExponentialCorrection");
alg->setProperty("C0", getPropertyValue("C0"));
alg->setProperty("C1", getPropertyValue("C1"));
} else if (correctionAlgorithm == "AutoDetect") {
// Figure out what to do from the instrument
try {
const auto corrVec = instrument->getStringParameter("correction");
if (corrVec.empty()) {
throw std::runtime_error(
"Could not find parameter 'correction' in "
"parameter file. Cannot auto detect the type of "
"correction.");
}
const std::string correctionStr = corrVec[0];
if (correctionStr == "polynomial") {
const auto polyVec = instrument->getStringParameter("polystring");
if (polyVec.empty())
throw std::runtime_error("Could not find parameter 'polystring' in "
"parameter file. Cannot apply polynomial "
"correction.");
alg->setProperty("CorrectionAlgorithm", "PolynomialCorrection");
alg->setProperty("Polynomial", polyVec[0]);
} else if (correctionStr == "exponential") {
const auto c0Vec = instrument->getStringParameter("C0");
if (c0Vec.empty())
throw std::runtime_error(
"Could not find parameter 'C0' in parameter "
"file. Cannot apply exponential correction.");
const auto c1Vec = instrument->getStringParameter("C1");
if (c1Vec.empty())
throw std::runtime_error(
"Could not find parameter 'C1' in parameter "
"file. Cannot apply exponential correction.");
alg->setProperty("C0", c0Vec[0]);
alg->setProperty("C1", c1Vec[0]);
}
alg->setProperty("NormalizeByIntegratedMonitors", false);
} catch (std::runtime_error &e) {
g_log.error() << e.what()
<< ". Polynomial correction will not be performed.";
alg->setProperty("CorrectionAlgorithm", "None");
}
} else {
alg->setProperty("CorrectionAlgorithm", "None");
}
}
/** Rebin and scale a workspace in Q.
*
* @param inputWS :: the workspace in Q
* @param theta :: the angle of this run
* @param params :: [output] rebin parameters
* @return :: the output workspace
*/
MatrixWorkspace_sptr
ReflectometryReductionOneAuto2::rebinAndScale(MatrixWorkspace_sptr inputWS,
const double theta,
std::vector<double> ¶ms) {
Property *qStepProp = getProperty("MomentumTransferStep");
double qstep;
if (!qStepProp->isDefault()) {
qstep = getProperty("MomentumTransferStep");
qstep = -qstep;
} else {
if (theta == 0.0) {
throw std::runtime_error(
"Theta determined from the detector positions is "
"0.0. Please provide a value for theta manually "
"or correct the detector position before running "
"this algorithm.");
}
IAlgorithm_sptr calcRes = createChildAlgorithm("NRCalculateSlitResolution");
calcRes->setProperty("Workspace", inputWS);
calcRes->setProperty("TwoTheta", 2 * theta);
calcRes->execute();
if (!calcRes->isExecuted()) {
g_log.error(
"NRCalculateSlitResolution failed. Workspace in Q will not be "
"rebinned. Please provide dQ/Q.");
return inputWS;
}
qstep = calcRes->getProperty("Resolution");
qstep = -qstep;
}
Property *qMin = getProperty("MomentumTransferMin");
Property *qMax = getProperty("MomentumTransferMax");
if (!qMin->isDefault() && !qMax->isDefault()) {
double qmin = getProperty("MomentumTransferMin");
double qmax = getProperty("MomentumTransferMax");
params.push_back(qmin);
params.push_back(qstep);
params.push_back(qmax);
} else {
params.push_back(qstep);
}
// Rebin
IAlgorithm_sptr algRebin = createChildAlgorithm("Rebin");
algRebin->initialize();
algRebin->setProperty("InputWorkspace", inputWS);
algRebin->setProperty("OutputWorkspace", inputWS);
algRebin->setProperty("Params", params);
algRebin->execute();
MatrixWorkspace_sptr IvsQ = algRebin->getProperty("OutputWorkspace");
// Scale (optional)
Property *scaleProp = getProperty("ScaleFactor");
if (!scaleProp->isDefault()) {
double scaleFactor = getProperty("ScaleFactor");
IAlgorithm_sptr algScale = createChildAlgorithm("Scale");
algScale->initialize();
algScale->setProperty("InputWorkspace", IvsQ);
algScale->setProperty("OutputWorkspace", IvsQ);
algScale->setProperty("Factor", 1.0 / scaleFactor);
algScale->execute();
IvsQ = algScale->getProperty("OutputWorkspace");
}
return IvsQ;
}
/** Check if input workspace is a group
*/
bool ReflectometryReductionOneAuto2::checkGroups() {
const std::string wsName = getPropertyValue("InputWorkspace");
try {
auto ws =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(wsName);
if (ws)
return true;
} catch (...) {
}
return false;
}
/** Process groups. Groups are processed differently depending on transmission
* runs and polarization analysis. If transmission run is a matrix workspace, it
* will be applied to each of the members in the input workspace group. If
* transmission run is a workspace group, the behaviour is different depending
* on polarization analysis. If polarization analysis is off (i.e.
* 'PolarizationAnalysis' is set to 'None') each item in the transmission group
* is associated with the corresponding item in the input workspace group. If
* polarization analysis is on (i.e. 'PolarizationAnalysis' is 'PA' or 'PNR')
* items in the transmission group will be summed to produce a matrix workspace
* that will be applied to each of the items in the input workspace group. See
* documentation of this algorithm for more details.
*/
bool ReflectometryReductionOneAuto2::processGroups() {
// this algorithm effectively behaves as MultiPeriodGroupAlgorithm
m_usingBaseProcessGroups = true;
// Get our input workspace group
auto group = AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(
getPropertyValue("InputWorkspace"));
// Get name of IvsQ workspace (native binning)
const std::string outputIvsQ = getPropertyValue("OutputWorkspace");
// Get name of IvsQ (native binning) workspace
const std::string outputIvsQBinned =
getPropertyValue("OutputWorkspaceBinned");
// Get name of IvsLam workspace
const std::string outputIvsLam =
getPropertyValue("OutputWorkspaceWavelength");
// Create a copy of ourselves
Algorithm_sptr alg =
createChildAlgorithm(name(), -1, -1, isLogging(), version());
alg->setChild(false);
alg->setRethrows(true);
// Copy all the non-workspace properties over
const std::vector<Property *> props = getProperties();
for (auto &prop : props) {
if (prop) {
IWorkspaceProperty *wsProp = dynamic_cast<IWorkspaceProperty *>(prop);
if (!wsProp)
alg->setPropertyValue(prop->name(), prop->value());
}
}
const bool polarizationAnalysisOn =
getPropertyValue("PolarizationAnalysis") != "None";
// Check if the transmission runs are groups or not
const std::string firstTrans = getPropertyValue("FirstTransmissionRun");
WorkspaceGroup_sptr firstTransG;
MatrixWorkspace_sptr firstTransSum;
if (!firstTrans.empty()) {
auto firstTransWS =
AnalysisDataService::Instance().retrieveWS<Workspace>(firstTrans);
firstTransG = boost::dynamic_pointer_cast<WorkspaceGroup>(firstTransWS);
if (!firstTransG) {
alg->setProperty("FirstTransmissionRun", firstTrans);
} else if (polarizationAnalysisOn) {
firstTransSum = sumTransmissionWorkspaces(firstTransG);
}
}
const std::string secondTrans = getPropertyValue("SecondTransmissionRun");
WorkspaceGroup_sptr secondTransG;
MatrixWorkspace_sptr secondTransSum;
if (!secondTrans.empty()) {
auto secondTransWS =
AnalysisDataService::Instance().retrieveWS<Workspace>(secondTrans);
secondTransG = boost::dynamic_pointer_cast<WorkspaceGroup>(secondTransWS);
if (!secondTransG) {
alg->setProperty("SecondTransmissionRun", secondTrans);
} else if (polarizationAnalysisOn) {
secondTransSum = sumTransmissionWorkspaces(secondTransG);
}
}
std::vector<std::string> IvsQGroup, IvsQUnbinnedGroup, IvsLamGroup;
// Execute algorithm over each group member
for (size_t i = 0; i < group->size(); ++i) {
const std::string IvsQName = outputIvsQ + "_" + std::to_string(i + 1);
const std::string IvsQBinnedName =
outputIvsQBinned + "_" + std::to_string(i + 1);
const std::string IvsLamName = outputIvsLam + "_" + std::to_string(i + 1);
if (firstTransG) {
if (!polarizationAnalysisOn)
alg->setProperty("FirstTransmissionRun",
firstTransG->getItem(i)->getName());
else
alg->setProperty("FirstTransmissionRun", firstTransSum);
}
if (secondTransG) {
if (!polarizationAnalysisOn)
alg->setProperty("SecondTransmissionRun",
secondTransG->getItem(i)->getName());
else
alg->setProperty("SecondTransmissionRun", secondTransSum);
}
alg->setProperty("InputWorkspace", group->getItem(i)->getName());
alg->setProperty("OutputWorkspace", IvsQName);
alg->setProperty("OutputWorkspaceBinned", IvsQBinnedName);
alg->setProperty("OutputWorkspaceWavelength", IvsLamName);
alg->execute();
IvsQGroup.push_back(IvsQName);
IvsQUnbinnedGroup.push_back(IvsQBinnedName);
IvsLamGroup.push_back(IvsLamName);
}
// Group the IvsQ and IvsLam workspaces
Algorithm_sptr groupAlg = createChildAlgorithm("GroupWorkspaces");
groupAlg->setChild(false);
groupAlg->setRethrows(true);
groupAlg->setProperty("InputWorkspaces", IvsLamGroup);
groupAlg->setProperty("OutputWorkspace", outputIvsLam);
groupAlg->execute();
groupAlg->setProperty("InputWorkspaces", IvsQGroup);
groupAlg->setProperty("OutputWorkspace", outputIvsQ);
groupAlg->execute();
groupAlg->setProperty("InputWorkspaces", IvsQUnbinnedGroup);
groupAlg->setProperty("OutputWorkspace", outputIvsQBinned);
groupAlg->execute();
// Set other properties so they can be updated in the Reflectometry interface
setPropertyValue("ThetaIn", alg->getPropertyValue("ThetaIn"));
setPropertyValue("MomentumTransferMin",
alg->getPropertyValue("MomentumTransferMin"));
setPropertyValue("MomentumTransferMax",
alg->getPropertyValue("MomentumTransferMax"));
setPropertyValue("MomentumTransferStep",
alg->getPropertyValue("MomentumTransferStep"));
setPropertyValue("ScaleFactor", alg->getPropertyValue("ScaleFactor"));
if (!polarizationAnalysisOn) {
// No polarization analysis. Reduction stops here
setPropertyValue("OutputWorkspace", outputIvsQ);
setPropertyValue("OutputWorkspaceBinned", outputIvsQBinned);
setPropertyValue("OutputWorkspaceWavelength", outputIvsLam);
return true;
}
applyPolarizationCorrection(outputIvsLam);
// Now we've overwritten the IvsLam workspaces, we'll need to recalculate
// the IvsQ ones
alg->setProperty("FirstTransmissionRun", "");
alg->setProperty("SecondTransmissionRun", "");
alg->setProperty("CorrectionAlgorithm", "None");
alg->setProperty("ProcessingInstructions", "0");
for (size_t i = 0; i < group->size(); ++i) {
const std::string IvsQName = outputIvsQ + "_" + std::to_string(i + 1);
const std::string IvsQBinnedName =
outputIvsQBinned + "_" + std::to_string(i + 1);
const std::string IvsLamName = outputIvsLam + "_" + std::to_string(i + 1);
alg->setProperty("InputWorkspace", IvsLamName);
alg->setProperty("OutputWorkspace", IvsQName);
alg->setProperty("OutputWorkspaceBinned", IvsQBinnedName);
alg->setProperty("OutputWorkspaceWavelength", IvsLamName);
alg->execute();
}
setPropertyValue("OutputWorkspace", outputIvsQ);
setPropertyValue("OutputWorkspaceBinned", outputIvsQBinned);
setPropertyValue("OutputWorkspaceWavelength", outputIvsLam);
return true;
}
/** Construct a polarization efficiencies workspace based on values of input
* properties.
*/
std::tuple<API::MatrixWorkspace_sptr, std::string, std::string>
ReflectometryReductionOneAuto2::getPolarizationEfficiencies() {
auto groupIvsLam = AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(
getPropertyValue("OutputWorkspaceWavelength"));
std::string const paMethod = getPropertyValue("PolarizationAnalysis");
Workspace_sptr workspace = groupIvsLam->getItem(0);
MatrixWorkspace_sptr efficiencies;
std::string correctionMethod;
std::string correctionOption;
if (paMethod == "ParameterFile") {
auto effAlg = createChildAlgorithm("ExtractPolarizationEfficiencies");
effAlg->setProperty("InputWorkspace", workspace);
effAlg->execute();
efficiencies = effAlg->getProperty("OutputWorkspace");
correctionMethod = effAlg->getPropertyValue("CorrectionMethod");
correctionOption = effAlg->getPropertyValue("CorrectionOption");
} else {
auto effAlg = createChildAlgorithm("CreatePolarizationEfficiencies");
effAlg->setProperty("InputWorkspace", workspace);
if (!isDefault("Pp")) {
effAlg->setProperty("Pp", getPropertyValue("Pp"));
}
if (!isDefault("Rho")) {
effAlg->setProperty("Rho", getPropertyValue("Rho"));
}
if (!isDefault("Ap")) {
effAlg->setProperty("Ap", getPropertyValue("Ap"));
}
if (!isDefault("Alpha")) {
effAlg->setProperty("Alpha", getPropertyValue("Alpha"));
}
effAlg->execute();
efficiencies = effAlg->getProperty("OutputWorkspace");
correctionMethod = "Fredrikze";
correctionOption = paMethod;
}
return std::make_tuple(efficiencies, correctionMethod, correctionOption);
}
/**
* Apply a polarization correction to workspaces in lambda.
* @param outputIvsLam :: Name of a workspace group to apply the correction to.
*/
void ReflectometryReductionOneAuto2::applyPolarizationCorrection(
std::string const &outputIvsLam) {
MatrixWorkspace_sptr efficiencies;
std::string correctionMethod;
std::string correctionOption;
std::tie(efficiencies, correctionMethod, correctionOption) =
getPolarizationEfficiencies();
Algorithm_sptr polAlg = createChildAlgorithm("PolarizationEfficiencyCor");
polAlg->setChild(false);
polAlg->setRethrows(true);
polAlg->setProperty("OutputWorkspace", outputIvsLam);
polAlg->setProperty("Efficiencies", efficiencies);
polAlg->setProperty("CorrectionMethod", correctionMethod);
if (correctionMethod == "Fredrikze") {
polAlg->setProperty("InputWorkspaceGroup", outputIvsLam);
polAlg->setProperty("PolarizationAnalysis", correctionOption);
} else {
throw std::invalid_argument("Unsupported polarization correction method: " +
correctionMethod);
}
polAlg->execute();
}
/**
* Sum transmission workspaces that belong to a workspace group
* @param transGroup : The transmission group containing the transmission runs
* @return :: A workspace pointer containing the sum of transmission workspaces
*/
MatrixWorkspace_sptr ReflectometryReductionOneAuto2::sumTransmissionWorkspaces(
WorkspaceGroup_sptr &transGroup) {
const std::string transSum = "trans_sum";
Workspace_sptr sumWS = transGroup->getItem(0)->clone();
/// For this step to appear in the history of the output workspaces I need to
/// set child to false and work with the ADS
auto plusAlg = createChildAlgorithm("Plus");
plusAlg->setChild(false);
plusAlg->initialize();
for (size_t item = 1; item < transGroup->size(); item++) {
plusAlg->setProperty("LHSWorkspace", sumWS);
plusAlg->setProperty("RHSWorkspace", transGroup->getItem(item));
plusAlg->setProperty("OutputWorkspace", transSum);
plusAlg->execute();
sumWS = AnalysisDataService::Instance().retrieve(transSum);
}
MatrixWorkspace_sptr result =
boost::dynamic_pointer_cast<MatrixWorkspace>(sumWS);
AnalysisDataService::Instance().remove(transSum);
return result;
}
} // namespace Algorithms
} // namespace Mantid