CalculateUMatrix.cpp 4.74 KB
Newer Older
1
2
3
4
#include "MantidCrystal/CalculateUMatrix.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
5
#include "MantidKernel/BoundedValidator.h"
6
7
8
9
10
11
12
13
14
15
16
17

namespace Mantid
{
namespace Crystal
{

  // Register the algorithm into the AlgorithmFactory
  DECLARE_ALGORITHM(CalculateUMatrix)

  using namespace Mantid::Kernel;
  using namespace Mantid::API;
  using namespace Mantid::DataObjects;
18
  using Mantid::Geometry::OrientedLattice;
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43


  //----------------------------------------------------------------------------------------------
  /** Constructor
   */
  CalculateUMatrix::CalculateUMatrix()
  {
  }
    
  //----------------------------------------------------------------------------------------------
  /** Destructor
   */
  CalculateUMatrix::~CalculateUMatrix()
  {
  }
  

  //----------------------------------------------------------------------------------------------

  //----------------------------------------------------------------------------------------------
  /** Initialize the algorithm's properties.
   */
  void CalculateUMatrix::init()
  {
    this->declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","",Direction::InOut), "An input workspace.");
44
    boost::shared_ptr<BoundedValidator<double> > mustBePositive = boost::make_shared<BoundedValidator<double> >();
45
    mustBePositive->setLower(0.0);
46
    boost::shared_ptr<BoundedValidator<double> > reasonable_angle = boost::make_shared<BoundedValidator<double> >();
47
48
49
    reasonable_angle->setLower(5.0);
    reasonable_angle->setUpper(175.0);
    // put in negative values, so user is forced to input all parameters. no shortcuts :)
50
51
52
53
54
55
    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");
56
57
58
59
60
61
62
63
64
65
66
67
68
69
  }

  //----------------------------------------------------------------------------------------------
  /** Execute the algorithm.
   */
  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");
    OrientedLattice o(a,b,c,alpha,beta,gamma);
70
    Matrix<double> B=o.getB();
71

72
    PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
73
74
    if (!ws) throw std::runtime_error("Problems reading the peaks workspace");

75
76
77
    size_t nIndexedpeaks=0;
    bool found2nc=false;
    V3D old(0,0,0);
78
    Matrix<double> Hi(4,4),Si(4,4),HS(4,4),zero(4,4);
79
80
81
    for (int i=0;i<ws->getNumberPeaks();i++)
    {
      Peak p=ws->getPeaks()[i];
Nick Draper's avatar
Nick Draper committed
82
83
84
      double H=p.getH();
      double K=p.getK();
      double L=p.getL();
85
86
      if(H*H+K*K+L*L>0)
      {
87
88
89
90
91
92
93
94
95
96
97
98
        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;
          }
        }
99
        V3D Qhkl=B*V3D(H,K,L);
100
101
102
103
104
        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.;

105
106
107
108
109
        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.;
110
111

        HS+=(Hi*Si);
112
113
      }
    }
114
115
116
    //check if enough peaks are indexed or if HS is 0
    if ((nIndexedpeaks<2) || (found2nc==false)) throw std::invalid_argument("Less then two non-colinear peaks indexed");
    if (HS==zero) throw std::invalid_argument("Something really bad happened");
117
118
119

    Matrix<double> Eval;
    Matrix<double> Diag;
120
    HS.Diagonalise(Eval,Diag);
121
    Eval.sortEigen(Diag);
122
    Mantid::Kernel::Quat qR(Eval[0][0],Eval[1][0],Eval[2][0],Eval[3][0]);//the first column corresponds to the highest eigenvalue
123
    DblMatrix U(qR.getRotation());
124
    o.setU(U);
125
126
127
128
129
130
131
132
133

    ws->mutableSample().setOrientedLattice(new OrientedLattice(o));

  }



} // namespace Mantid
} // namespace Crystal