CalculateUMatrix.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCrystal/CalculateUMatrix.h"
#include "MantidAPI/Sample.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/BoundedValidator.h"
 
namespace Mantid::Crystal {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(CalculateUMatrix)
 
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using Mantid::Geometry::OrientedLattice;
 
void CalculateUMatrix::init() {
  this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
                        "An input workspace.");
  std::shared_ptr<BoundedValidator<double>> mustBePositive = std::make_shared<BoundedValidator<double>>();
  mustBePositive->setLower(0.0);
  std::shared_ptr<BoundedValidator<double>> reasonable_angle = std::make_shared<BoundedValidator<double>>();
  reasonable_angle->setLower(5.0);
  reasonable_angle->setUpper(175.0);
  // put in negative values, so user is forced to input all parameters. no
  // shortcuts :)
  this->declareProperty("a", -1.0, mustBePositive, "Lattice parameter a");
  this->declareProperty("b", -1.0, mustBePositive, "Lattice parameter b");
  this->declareProperty("c", -1.0, mustBePositive, "Lattice parameter c");
  this->declareProperty("alpha", -1.0, reasonable_angle, "Lattice parameter alpha");
  this->declareProperty("beta", -1.0, reasonable_angle, "Lattice parameter beta");
  this->declareProperty("gamma", -1.0, reasonable_angle, "Lattice parameter gamma");
}
 
void CalculateUMatrix::exec() {
  double a = this->getProperty("a");
  double b = this->getProperty("b");
  double c = this->getProperty("c");
  double alpha = this->getProperty("alpha");
  double beta = this->getProperty("beta");
  double gamma = this->getProperty("gamma");
  auto lattice = std::make_unique<OrientedLattice>(a, b, c, alpha, beta, gamma);
  Matrix<double> B = lattice->getB();
 
  IPeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
  if (!ws)
    throw std::runtime_error("Problems reading the peaks workspace");
 
  size_t nIndexedpeaks = 0;
  bool found2nc = false;
  V3D old(0, 0, 0);
  Matrix<double> Hi(4, 4), Si(4, 4), HS(4, 4), zero(4, 4);
  for (int i = 0; i < ws->getNumberPeaks(); i++) {
    Mantid::Geometry::IPeak &p = ws->getPeak(i);
    double H = p.getH();
    double K = p.getK();
    double L = p.getL();
    if (H * H + K * K + L * L > 0) {
      nIndexedpeaks++;
      if (!found2nc) {
        if (nIndexedpeaks == 1) {
          old = V3D(H, K, L);
        } else {
          if (!old.coLinear(V3D(0, 0, 0), V3D(H, K, L)))
            found2nc = true;
        }
      }
      V3D Qhkl = B * V3D(H, K, L);
      Hi[0][0] = 0.;
      Hi[0][1] = -Qhkl.X();
      Hi[0][2] = -Qhkl.Y();
      Hi[0][3] = -Qhkl.Z();
      Hi[1][0] = Qhkl.X();
      Hi[1][1] = 0.;
      Hi[1][2] = Qhkl.Z();
      Hi[1][3] = -Qhkl.Y();
      Hi[2][0] = Qhkl.Y();
      Hi[2][1] = -Qhkl.Z();
      Hi[2][2] = 0.;
      Hi[2][3] = Qhkl.X();
      Hi[3][0] = Qhkl.Z();
      Hi[3][1] = Qhkl.Y();
      Hi[3][2] = -Qhkl.X();
      Hi[3][3] = 0.;
 
      V3D Qgon = p.getQSampleFrame();
      Si[0][0] = 0.;
      Si[0][1] = -Qgon.X();
      Si[0][2] = -Qgon.Y();
      Si[0][3] = -Qgon.Z();
      Si[1][0] = Qgon.X();
      Si[1][1] = 0.;
      Si[1][2] = -Qgon.Z();
      Si[1][3] = Qgon.Y();
      Si[2][0] = Qgon.Y();
      Si[2][1] = Qgon.Z();
      Si[2][2] = 0.;
      Si[2][3] = -Qgon.X();
      Si[3][0] = Qgon.Z();
      Si[3][1] = -Qgon.Y();
      Si[3][2] = Qgon.X();
      Si[3][3] = 0.;
 
      HS += (Hi * Si);
    }
  }
  // check if enough peaks are indexed or if HS is 0
  if ((nIndexedpeaks < 2) || !found2nc)
    throw std::invalid_argument("Less then two non-colinear peaks indexed");
  if (HS == zero)
    throw std::invalid_argument("Something really bad happened");
 
  Matrix<double> Eval;
  Matrix<double> Diag;
  HS.Diagonalise(Eval, Diag);
  Eval.sortEigen(Diag);
  Mantid::Kernel::Quat qR(Eval[0][0], Eval[1][0], Eval[2][0],
                          Eval[3][0]); // the first column corresponds to the highest eigenvalue
  DblMatrix U(qR.getRotation());
  lattice->setU(U);
 
  ws->mutableSample().setOrientedLattice(std::move(lattice));
}
 
} // namespace Mantid::Crystal