CreatePolarizationEfficiencies.cpp

#include "MantidDataHandling/CreatePolarizationEfficiencies.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceHistory.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/Unit.h"
#include "MantidGeometry/Instrument.h"

#include <boost/shared_ptr.hpp>

#include <algorithm>

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::Geometry;

namespace {

double calculatePolynomial(std::vector<double> const &coefficients, double x) {
  double polynomial = coefficients[0];
  double xPow = 1.0;
  // Build up the polynomial in ascending powers of x
  for (size_t i = 1; i < coefficients.size(); ++i) {
    xPow *= x;
    polynomial += coefficients[i] * xPow;
  }
  return polynomial;
}

} // namespace

namespace Mantid {
namespace DataHandling {

DECLARE_ALGORITHM(CreatePolarizationEfficiencies)

const std::string CreatePolarizationEfficiencies::name() const {
  return "CreatePolarizationEfficiencies";
}

int CreatePolarizationEfficiencies::version() const { return 1; }

const std::string CreatePolarizationEfficiencies::summary() const {
  return "Converts polynomial factors to histograms with polarization "
         "efficiencies.";
}

const std::vector<std::string> CreatePolarizationEfficiencies::seeAlso() const {
  return {"JoinISISPolarizationEfficiencies",
          "LoadISISPolarizationEfficiencies", "PolarizationEfficiencyCor"};
}

void CreatePolarizationEfficiencies::init() {
  declareProperty(make_unique<WorkspaceProperty<Mantid::API::MatrixWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
                  "An input workspace to use the x-values from.");

  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.");

  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>(P1, Direction::Input),
      "Polarizer efficiency.");

  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>(P2, Direction::Input),
      "Analyzer efficiency.");

  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>(F1, Direction::Input),
      "Polarizer flipper efficiency.");

  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>(F2, Direction::Input),
      "Analyzer flipper efficiency.");

  initOutputWorkspace();
}

/// Create the efficiencies workspace given names of input properties.
/// @param labels :: Names of efficiencies which to include in the output
/// workspace.
MatrixWorkspace_sptr CreatePolarizationEfficiencies::createEfficiencies(
    std::vector<std::string> const &labels) {

  std::vector<std::vector<double>> polynomialCoefficients;

  for (auto const &label : labels) {
    polynomialCoefficients.emplace_back<std::vector<double>>(
        getProperty(label));
  }

  MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
  auto sharedInX = inWS->sharedX(0);

  MatrixWorkspace_sptr outWS = WorkspaceFactory::Instance().create(
      inWS, labels.size(), sharedInX->size(), inWS->blocksize());
  auto axis1 = new TextAxis(labels.size());
  outWS->replaceAxis(1, axis1);
  outWS->getAxis(0)->setUnit(inWS->getAxis(0)->unit()->unitID());

  auto const x = inWS->points(0);
  std::vector<double> y(x.size());
  for (size_t i = 0; i < labels.size(); ++i) {
    outWS->setSharedX(i, sharedInX);
    auto const &coefficients = polynomialCoefficients[i];
    std::transform(x.begin(), x.end(), y.begin(), [&coefficients](double v) {
      return calculatePolynomial(coefficients, v);
    });
    outWS->mutableY(i) = y;
    axis1->setLabel(i, labels[i]);
  }

  return outWS;
}

} // namespace Algorithms
} // namespace Mantid