IndexPeaks.cpp

#include "MantidCrystal/IndexPeaks.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidAPI/Sample.h"

namespace Mantid {
namespace Crystal {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IndexPeaks)

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

/** Initialize the algorithm's properties.
 */
void IndexPeaks::init() {
  this->declareProperty(make_unique<WorkspaceProperty<PeaksWorkspace> >(
                            "PeaksWorkspace", "", Direction::InOut),
                        "Input Peaks Workspace");

  auto mustBePositive = boost::make_shared<BoundedValidator<double> >();
  mustBePositive->setLower(0.0);

  this->declareProperty(
      make_unique<PropertyWithValue<double> >("Tolerance", 0.15, mustBePositive,
                                              Direction::Input),
      "Indexing Tolerance (0.15)");

  this->declareProperty(
      make_unique<PropertyWithValue<int> >("NumIndexed", 0, Direction::Output),
      "Gets set with the number of indexed peaks.");

  this->declareProperty(make_unique<PropertyWithValue<double> >(
                            "AverageError", 0.0, Direction::Output),
                        "Gets set with the average HKL indexing error.");

  this->declareProperty("RoundHKLs", true,
                        "Round H, K and L values to integers");

  this->declareProperty("CommonUBForAll", false,
                        "Index all orientations with a common UB");
}

/** Execute the algorithm.
 */
void IndexPeaks::exec() {
  PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
  if (!ws) {
    throw std::runtime_error("Could not read the peaks workspace");
  }

  OrientedLattice o_lattice = ws->mutableSample().getOrientedLattice();
  const Matrix<double> &UB = o_lattice.getUB();

  if (!IndexingUtils::CheckUB(UB)) {
    throw std::runtime_error(
        "ERROR: The stored UB is not a valid orientation matrix");
  }

  bool round_hkls = this->getProperty("RoundHKLs");
  bool commonUB = this->getProperty("CommonUBForAll");

  std::vector<Peak> &peaks = ws->getPeaks();
  size_t n_peaks = ws->getNumberPeaks();
  int total_indexed = 0;
  double average_error;
  double total_error = 0;
  double tolerance = this->getProperty("Tolerance");

  if (commonUB) {
    std::vector<V3D> miller_indices;
    std::vector<V3D> q_vectors;

    q_vectors.reserve(n_peaks);
    for (size_t i = 0; i < n_peaks; i++) {
      q_vectors.push_back(peaks[i].getQSampleFrame());
    }

    total_indexed = IndexingUtils::CalculateMillerIndices(
        UB, q_vectors, tolerance, miller_indices, average_error);

    for (size_t i = 0; i < n_peaks; i++) {
      peaks[i].setHKL(miller_indices[i]);
    }
  } else {
    // get list of run numbers in this peaks workspace
    std::vector<int> run_numbers;
    for (size_t i = 0; i < n_peaks; i++) {
      int run = peaks[i].getRunNumber();
      bool found = false;
      size_t k = 0;
      while (k < run_numbers.size() && !found) {
        if (run == run_numbers[k])
          found = true;
        else
          k++;
      }
      if (!found)
        run_numbers.push_back(run);
    }

    // index the peaks for each run separately, using a UB matrix optimized for
    // that run

    for (size_t run_index = 0; run_index < run_numbers.size(); run_index++) {
      std::vector<V3D> miller_indices;
      std::vector<V3D> q_vectors;

      int run = run_numbers[run_index];
      for (size_t i = 0; i < n_peaks; i++) {
        if (peaks[i].getRunNumber() == run)
          q_vectors.push_back(peaks[i].getQSampleFrame());
      }

      Matrix<double> tempUB(UB);

      int num_indexed = 0;
      int original_indexed = 0;
      double original_error = 0;
      original_indexed = IndexingUtils::CalculateMillerIndices(
          UB, q_vectors, tolerance, miller_indices, original_error);

      IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded for
                                                // Optimize_UB to work
      num_indexed = original_indexed;
      average_error = original_error;

      bool done = false;
      if (num_indexed < 3) // can't optimize without at least 3
      {                    // peaks
        done = true;
      }

      int iteration = 0;
      while (iteration < 4 && !done) // try repeatedly optimizing 4 times
      {                              // which is usually sufficient
        try {
          IndexingUtils::Optimize_UB(tempUB, miller_indices, q_vectors);
        }
        catch (...) // If there is any problem, such as too few
        {           // independent peaks, just use the original UB
          tempUB = UB;
          done = true;
        }

        num_indexed = IndexingUtils::CalculateMillerIndices(
            UB, q_vectors, tolerance, miller_indices, average_error);

        IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded for
                                                  // Optimize_UB to work

        if (num_indexed < original_indexed) // just use the original UB
        {
          num_indexed = original_indexed;
          average_error = original_error;
          done = true;
        }

        iteration++;
      }

      if (!round_hkls) // If user wants fractional hkls, recalculate them
      {
        num_indexed = IndexingUtils::CalculateMillerIndices(
            UB, q_vectors, tolerance, miller_indices, average_error);
      }

      total_indexed += num_indexed;
      total_error += average_error * num_indexed;

      // tell the user how many were indexed in each run
      if (run_numbers.size() > 1) {
        g_log.notice() << "Run " << run << ": indexed " << num_indexed
                       << " Peaks out of " << q_vectors.size()
                       << " with tolerance of " << tolerance << '\n';
        g_log.notice() << "Average error in h,k,l for indexed peaks =  "
                       << average_error << '\n';
      }

      size_t miller_index_counter = 0;
      for (size_t i = 0; i < n_peaks; i++) {
        if (peaks[i].getRunNumber() == run) {
          peaks[i].setHKL(miller_indices[miller_index_counter]);
          miller_index_counter++;
        }
      }
    }

    if (total_indexed > 0)
      average_error = total_error / total_indexed;
    else
      average_error = 0;
  }

  // tell the user how many were indexed overall and the overall average error
  g_log.notice() << "ALL Runs: indexed " << total_indexed << " Peaks out of "
                 << n_peaks << " with tolerance of " << tolerance << '\n';
  g_log.notice() << "Average error in h,k,l for indexed peaks =  "
                 << average_error << '\n';

  // Save output properties
  this->setProperty("NumIndexed", total_indexed);
  this->setProperty("AverageError", average_error);
  // Show the lattice parameters
  g_log.notice() << o_lattice << "\n";
}

} // namespace Mantid
} // namespace Crystal