#include "MantidMDAlgorithms/ConvertToReflectometryQ.h" #include "MantidAPI/IEventWorkspace.h" #include "MantidAPI/ITableWorkspace.h" #include "MantidAPI/WorkspaceValidators.h" #include "MantidAPI/Progress.h" #include "MantidDataObjects/EventWorkspace.h" #include "MantidDataObjects/TableWorkspace.h" #include "MantidDataObjects/Workspace2D.h" #include "MantidKernel/ArrayProperty.h" #include "MantidKernel/EnabledWhenProperty.h" #include "MantidKernel/ListValidator.h" #include "MantidKernel/TimeSeriesProperty.h" #include "MantidDataObjects/MDEventWorkspace.h" #include "MantidDataObjects/MDEventFactory.h" #include "MantidMDAlgorithms/ReflectometryTransformKiKf.h" #include "MantidMDAlgorithms/ReflectometryTransformQxQz.h" #include "MantidMDAlgorithms/ReflectometryTransformP.h" #include <boost/optional.hpp> #include <boost/shared_ptr.hpp> #include <boost/make_shared.hpp> using namespace Mantid::API; using namespace Mantid::Kernel; using namespace Mantid::Geometry; using namespace Mantid::DataObjects; /*Non member helpers*/ namespace { /* Transform to q-space label: @return: associated id/label */ std::string qSpaceTransform() { return "Q (lab frame)"; } /* Transform to p-space label: @return: associated id/label */ std::string pSpaceTransform() { return "P (lab frame)"; } /* Transform to k-space label: @return: associated id/label */ std::string kSpaceTransform() { return "K (incident, final)"; } /* Center point rebinning @return: associated id/label */ std::string centerTransform() { return "Centre"; } /* Normalised polygon rebinning @return: associated id/label */ std::string normPolyTransform() { return "NormalisedPolygon"; } /* Check that the input workspace is of the correct type. @param: inputWS: The input workspace. @throw: runtime_error if the units do not appear to be correct/compatible with the algorithm. */ void checkInputWorkspace(Mantid::API::MatrixWorkspace_const_sptr inputWs) { auto spectraAxis = inputWs->getAxis(1); const std::string label = spectraAxis->unit()->label(); const std::string expectedLabel = "degrees"; if (expectedLabel != label) { std::string message = "Spectra axis should have units of degrees. Instead found: " + label; throw std::runtime_error(message); } } /* Check the extents. @param extents: A vector containing all the extents. @throw: runtime_error if the extents appear to be incorrect. */ void checkExtents(const std::vector<double> &extents) { if (extents.size() != 4) { throw std::runtime_error( "Four comma separated extents inputs should be provided"); } if ((extents[0] >= extents[1]) || (extents[2] >= extents[3])) { throw std::runtime_error( "Extents must be provided min, max with min less than max!"); } } /* Check the incident theta inputs. @param bUseOwnIncidentTheta: True if the user has requested to provide their own incident theta value. @param theta: The proposed incident theta value provided by the user. @throw: logic_error if the theta value is out of range. */ void checkCustomThetaInputs(const bool bUseOwnIncidentTheta, const double &theta) { if (bUseOwnIncidentTheta) { if (theta < 0 || theta > 90) { throw std::logic_error("Overriding incident theta is out of range"); } } } /* General check for the indient theta. @param theta: The proposed incident theta value. @throw: logic_error if the theta value is out of range. */ void checkIncidentTheta(const double &theta) { if (theta < 0 || theta > 90) { throw std::logic_error("Incident theta is out of range"); } } /* Check for the output dimensionality. @param outputDimensions : requested output dimensionality @throw: runtime_errror if the dimensionality is not supported. */ void checkOutputDimensionalityChoice(const std::string &outputDimensions) { if (outputDimensions != qSpaceTransform() && outputDimensions != kSpaceTransform() && outputDimensions != pSpaceTransform()) { throw std::runtime_error("Unknown transformation"); } } } namespace Mantid { namespace MDAlgorithms { // Register the algorithm into the AlgorithmFactory DECLARE_ALGORITHM(ConvertToReflectometryQ) //---------------------------------------------------------------------------------------------- /** Constructor */ ConvertToReflectometryQ::ConvertToReflectometryQ() {} //---------------------------------------------------------------------------------------------- /** Destructor */ ConvertToReflectometryQ::~ConvertToReflectometryQ() {} //---------------------------------------------------------------------------------------------- /// Algorithm's name for identification. @see Algorithm::name const std::string ConvertToReflectometryQ::name() const { return "ConvertToReflectometryQ"; } /// Algorithm's version for identification. @see Algorithm::version int ConvertToReflectometryQ::version() const { return 1; } /// Algorithm's category for identification. @see Algorithm::category const std::string ConvertToReflectometryQ::category() const { return "Reflectometry"; } //---------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------- /** Initialize the algorithm's properties. */ void ConvertToReflectometryQ::init() { auto compositeValidator = boost::make_shared<CompositeValidator>(); compositeValidator->add( boost::make_shared<API::WorkspaceUnitValidator>("Wavelength")); compositeValidator->add(boost::make_shared<API::HistogramValidator>()); declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "", Direction::Input, compositeValidator), "An input workspace in wavelength"); std::vector<std::string> propOptions; propOptions.push_back(qSpaceTransform()); propOptions.push_back(pSpaceTransform()); propOptions.push_back(kSpaceTransform()); declareProperty( "OutputDimensions", "Q (lab frame)", boost::make_shared<StringListValidator>(propOptions), "What will be the dimensions of the output workspace?\n" " Q (lab frame): Wave-vector change of the lattice in the lab frame.\n" " P (lab frame): Momentum in the sample frame.\n" " K initial and final vectors in the z plane."); std::vector<std::string> transOptions; transOptions.push_back(centerTransform()); transOptions.push_back(normPolyTransform()); declareProperty("Method", centerTransform(), boost::make_shared<StringListValidator>(transOptions), "What method should be used for the axis transformation?\n" " Centre: Use center point rebinning.\n" " NormalisedPolygon: Use normalised polygon rebinning."); declareProperty( new Kernel::PropertyWithValue<bool>("OverrideIncidentTheta", false), "Use the provided incident theta value."); declareProperty( new PropertyWithValue<double>("IncidentTheta", -1), "An input incident theta value specified in degrees." "Optional input value for the incident theta specified in degrees."); std::vector<double> extents(4, 0); extents[0] = -50; extents[1] = +50; extents[2] = -50; extents[3] = +50; declareProperty(new ArrayProperty<double>("Extents", extents), "A comma separated list of min, max for each dimension. " "Takes four values in the form dim_0_min, dim_0_max, " "dim_1_min, dim_1_max,\n" "specifying the extents of each dimension. Optional, default " "+-50 in each dimension."); setPropertySettings("IncidentTheta", new Kernel::EnabledWhenProperty("OverrideIncidentTheta", IS_EQUAL_TO, "1")); declareProperty( new Kernel::PropertyWithValue<bool>("OutputAsMDWorkspace", true), "Generate the output as a MDWorkspace, otherwise a Workspace2D is " "returned."); declareProperty(new WorkspaceProperty<IMDWorkspace>("OutputWorkspace", "", Direction::Output), "Output 2D Workspace."); declareProperty(new WorkspaceProperty<ITableWorkspace>("OutputVertexes", "", Direction::Output), "Output TableWorkspace with vertex information. See " "DumpVertexes property."); declareProperty(new Kernel::PropertyWithValue<int>("NumberBinsQx", 100), "The number of bins along the qx axis. Optional and only " "applies to 2D workspaces. Defaults to 100."); declareProperty(new Kernel::PropertyWithValue<int>("NumberBinsQz", 100), "The number of bins along the qx axis. Optional and only " "applies to 2D workspaces. Defaults to 100."); declareProperty(new Kernel::PropertyWithValue<bool>("DumpVertexes", false), "If set, with 2D rebinning, the intermediate vertexes for " "each polygon will be written out for debugging purposes. " "Creates a second output table workspace."); setPropertySettings( "NumberBinsQx", new EnabledWhenProperty("OutputAsMDWorkspace", IS_NOT_DEFAULT)); setPropertySettings( "NumberBinsQz", new EnabledWhenProperty("OutputAsMDWorkspace", IS_NOT_DEFAULT)); // Create box controller properties. this->initBoxControllerProps("2,2", 50, 10); // Only show box controller properties when a md workspace is returned. setPropertySettings( "SplitInto", new EnabledWhenProperty("OutputAsMDWorkspace", IS_DEFAULT)); setPropertySettings("SplitThreshold", new EnabledWhenProperty( "OutputAsMDWorkspace", IS_DEFAULT)); setPropertySettings( "MaxRecursionDepth", new EnabledWhenProperty("OutputAsMDWorkspace", IS_DEFAULT)); } //---------------------------------------------------------------------------------------------- /** Execute the algorithm. */ void ConvertToReflectometryQ::exec() { Mantid::API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace"); const bool bUseOwnIncidentTheta = getProperty("OverrideIncidentTheta"); const std::vector<double> extents = getProperty("Extents"); double incidentTheta = getProperty("IncidentTheta"); const std::string outputDimensions = getPropertyValue("OutputDimensions"); const std::string transMethod = getPropertyValue("Method"); const bool outputAsMDWorkspace = getProperty("OutputAsMDWorkspace"); const int numberOfBinsQx = getProperty("NumberBinsQx"); const int numberOfBinsQz = getProperty("NumberBinsQz"); // Validation of input parameters checkInputWorkspace(inputWs); checkExtents(extents); checkCustomThetaInputs(bUseOwnIncidentTheta, incidentTheta); checkOutputDimensionalityChoice(outputDimensions); // TODO: This check can be // retired as soon as all // transforms have been // Extract the incient theta angle from the logs if a user provided one is not // given. if (!bUseOwnIncidentTheta) { const Mantid::API::Run &run = inputWs->run(); try { Property *p = run.getLogData("stheta"); auto incidentThetas = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<double> *>(p); if (!incidentThetas) { throw std::runtime_error("stheta log not found"); } incidentTheta = incidentThetas->valuesAsVector().back(); // Not quite sure what to do // with the time series for // stheta checkIncidentTheta(incidentTheta); std::stringstream stream; stream << "Extracted initial theta value of: " << incidentTheta; g_log.information(stream.str()); } catch (Exception::NotFoundError &) { throw std::runtime_error( "The input workspace does not have a stheta log value."); } } // Min max extent values. const double dim0min = extents[0]; const double dim0max = extents[1]; const double dim1min = extents[2]; const double dim1max = extents[3]; BoxController_sptr bc = boost::make_shared<BoxController>(2); this->setBoxController(bc); // Select the transform strategy. ReflectometryTransform_sptr transform; if (outputDimensions == qSpaceTransform()) { transform = boost::make_shared<ReflectometryTransformQxQz>( dim0min, dim0max, dim1min, dim1max, incidentTheta, numberOfBinsQx, numberOfBinsQz); } else if (outputDimensions == pSpaceTransform()) { transform = boost::make_shared<ReflectometryTransformP>( dim0min, dim0max, dim1min, dim1max, incidentTheta, numberOfBinsQx, numberOfBinsQz); } else { transform = boost::make_shared<ReflectometryTransformKiKf>( dim0min, dim0max, dim1min, dim1max, incidentTheta, numberOfBinsQx, numberOfBinsQz); } IMDWorkspace_sptr outputWS; TableWorkspace_sptr vertexes = boost::make_shared<Mantid::DataObjects::TableWorkspace>(); Progress transSelectionProg(this, 0.0, 0.1, 2); if (outputAsMDWorkspace) { transSelectionProg.report("Choosing Transformation"); if (transMethod == centerTransform()) { auto outputMDWS = transform->executeMD(inputWs, bc); Progress transPerformProg(this, 0.1, 0.7, 5); transPerformProg.report("Performed transformation"); // Copy ExperimentInfo (instrument, run, sample) to the output WS ExperimentInfo_sptr ei(inputWs->cloneExperimentInfo()); outputMDWS->addExperimentInfo(ei); outputWS = outputMDWS; } else if (transMethod == normPolyTransform()) { Progress transPerformProg(this, 0.1, 0.7, 5); const bool dumpVertexes = this->getProperty("DumpVertexes"); auto vertexesTable = vertexes; // perform the normalised polygon transformation transPerformProg.report("Performing Transformation"); auto normPolyTrans = transform->executeNormPoly( inputWs, vertexesTable, dumpVertexes, outputDimensions); // copy any experiment info from input workspace normPolyTrans->copyExperimentInfoFrom(inputWs.get()); // produce MDHistoWorkspace from normPolyTrans workspace. Progress outputToMDProg(this, 0.7, 0.75, 10); auto outputMDWS = transform->executeMDNormPoly(normPolyTrans); ExperimentInfo_sptr ei(normPolyTrans->cloneExperimentInfo()); outputMDWS->addExperimentInfo(ei); outputWS = outputMDWS; outputToMDProg.report("Successfully output to MD"); } else { throw std::runtime_error("Unknown rebinning method: " + transMethod); } } else if (transMethod == normPolyTransform()) { transSelectionProg.report("Choosing Transformation"); Progress transPerformProg(this, 0.1, 0.7, 5); const bool dumpVertexes = this->getProperty("DumpVertexes"); auto vertexesTable = vertexes; // perform the normalised polygon transformation transPerformProg.report("Performing Transformation"); auto output2DWS = transform->executeNormPoly( inputWs, vertexesTable, dumpVertexes, outputDimensions); // copy any experiment info from input workspace output2DWS->copyExperimentInfoFrom(inputWs.get()); outputWS = output2DWS; transPerformProg.report("Transformation Complete"); } else if (transMethod == centerTransform()) { transSelectionProg.report("Choosing Transformation"); Progress transPerformProg(this, 0.1, 0.7, 5); transPerformProg.report("Performing Transformation"); auto output2DWS = transform->execute(inputWs); output2DWS->copyExperimentInfoFrom(inputWs.get()); outputWS = output2DWS; } else { throw std::runtime_error("Unknown rebinning method: " + transMethod); } // Execute the transform and bind to the output. setProperty("OutputWorkspace", outputWS); setProperty("OutputVertexes", vertexes); Progress setPropertyProg(this, 0.8, 1.0, 2); setPropertyProg.report("Success"); } } // namespace Mantid } // namespace MDAlgorithms