SelectCellWithForm.cpp 7.03 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "MantidCrystal/SelectCellWithForm.h"
#include "MantidCrystal/IndexPeaks.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/ScalarUtils.h"
#include "MantidGeometry/Crystal/ReducedCell.h"
#include "MantidGeometry/Crystal/ConventionalCell.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
10
#include "MantidKernel/BoundedValidator.h"
11
12
13
14
15
16
17
18
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
44
45
46
47
48
#include <boost/lexical_cast.hpp>
#include <cstdio>

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

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

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


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

49
    auto mustBePositive = boost::make_shared<BoundedValidator<int> >();
50
51
52
    mustBePositive->setLower(1);

    this->declareProperty(
53
          new PropertyWithValue<int>("FormNumber",0,mustBePositive,Direction::Input),
54
55
56
57
58
59
60
61
62
         "Form number for the desired cell");
    this->declareProperty( "Apply", false, "Update UB and re-index the peaks");
    this->declareProperty( "Tolerance", 0.12, "Indexing Tolerance");

    this->declareProperty(new PropertyWithValue<int>( "NumIndexed", 0,
          Direction::Output), "The number of indexed peaks if apply==true.");

    this->declareProperty(new PropertyWithValue<double>( "AverageError", 0.0,
          Direction::Output), "The average HKL indexing error if apply==true.");
63

64
    this->declareProperty( "AllowPermutations", true,
65
                            "Allow permutations of conventional cells" );
66
67
  }

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
  Kernel::Matrix<double> SelectCellWithForm::DetermineErrors( std::vector<double> &sigabc, const Kernel::Matrix<double> &UB,
                                             const PeaksWorkspace_sptr &ws, double tolerance)
  {

        std::vector<V3D> miller_ind;
        std::vector<V3D> q_vectors;
        std::vector<V3D>q_vectors0;
        int npeaks = ws->getNumberPeaks();
        double fit_error;
        miller_ind.reserve( npeaks );
        q_vectors.reserve( npeaks );
        q_vectors0.reserve(npeaks);
        for( int i=0;i<npeaks;i++)
          q_vectors0.push_back(ws->getPeak(i).getQSampleFrame());

        Kernel::Matrix<double>newUB1(3,3);
        IndexingUtils::GetIndexedPeaks(UB, q_vectors0, tolerance,
                            miller_ind, q_vectors, fit_error );
        IndexingUtils::Optimize_UB(newUB1, miller_ind,q_vectors,sigabc);

        int nindexed_old = (int)q_vectors.size();
        int nindexed_new = IndexingUtils::NumberIndexed( newUB1,q_vectors0,tolerance);
        bool latErrorsValid =true;
        if( nindexed_old < .8*nindexed_new || .8*nindexed_old >  nindexed_new )
           latErrorsValid = false;
        else
        {
          double maxDiff=0;
          double maxEntry =0;
          for( int row=0; row <3; row++)
            for( int col=0; col< 3; col++)
            {
              double diff= fabs( UB[row][col]- newUB1[row][col]);
              double V = std::max<double>(fabs(UB[row][col]), fabs(newUB1[row][col]));
              if( diff > maxDiff)
                maxDiff = diff;
              if(V > maxEntry)
                maxEntry = V;
            }
          if( maxEntry==0 || maxDiff/maxEntry >.1)
            latErrorsValid=false;
        }

        if( !latErrorsValid)
        {
          for( size_t i = 0; i < sigabc.size(); i++ )
            sigabc[i]=0;
          return UB;

        }else

          return newUB1;

  }

123
124
125
126
127
  //--------------------------------------------------------------------------
  /** Execute the algorithm.
   */
  void SelectCellWithForm::exec()
  {
128
    PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
129
130
131
132
133
134
135
136
    if (!ws) 
    { 
      throw std::runtime_error("Could not read the peaks workspace");
    }

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

137
138
    bool   allowPermutations        = this->getProperty("AllowPermutations");

139
140
141
142
143
144
145
146
147
148
    if ( ! IndexingUtils::CheckUB( UB ) )
    {
       throw std::runtime_error(
             "ERROR: The stored UB is not a valid orientation matrix");
    }

    int    form_num  = this->getProperty("FormNumber");
    bool   apply     = this->getProperty("Apply");
    double tolerance = this->getProperty("Tolerance");

149
    ConventionalCell info = ScalarUtils::GetCellForForm( UB, form_num, allowPermutations );
150
151
152

    DblMatrix newUB = info.GetNewUB();

153
154
    std::string message = info.GetDescription() + " Lat Par:" +
                          IndexingUtils::GetLatticeParameterString( newUB );
155
156
157

    g_log.notice( std::string(message) );

158
159
160
161
162
    Kernel::Matrix<double>T(UB);
    T.Invert();
    T = newUB * T;
    g_log.notice() << "Transformation Matrix =  " << T.str() << std::endl;

163
    if ( apply )
164
    {
165
166
167
168
169
    //----------------------------------- Try to optimize(LSQ) to find lattice errors ------------------------
    //                       UB matrix may NOT have been found by unconstrained least squares optimization

     //----------------------------------------------
    o_lattice.setUB( newUB );
170
171
172
173
    std::vector<double> sigabc(6);
    DetermineErrors(sigabc,newUB,ws, tolerance);

    o_lattice.setError( sigabc[0],sigabc[1],sigabc[2],sigabc[3],sigabc[4],sigabc[5]);
174

175
    ws->mutableSample().setOrientedLattice( new OrientedLattice(o_lattice) );
176

177
178
    std::vector<Peak> &peaks = ws->getPeaks();
    size_t n_peaks = ws->getNumberPeaks();
179

180
181
182
183
184
185
186
187
    int    num_indexed   = 0;
    double average_error = 0.0;
    std::vector<V3D> miller_indices;
    std::vector<V3D> q_vectors;
    for ( size_t i = 0; i < n_peaks; i++ )
    {
      q_vectors.push_back( peaks[i].getQSampleFrame() );
    }
188

189
190
191
192
    num_indexed = IndexingUtils::CalculateMillerIndices( newUB, q_vectors,
                                                         tolerance,
                                                         miller_indices,
                                                         average_error );
193

194
195
196
197
    for ( size_t i = 0; i < n_peaks; i++ )
    {
      peaks[i].setHKL( miller_indices[i] );
    }
198

199
200
201
    // Tell the user what happened.
    g_log.notice() << "Re-indexed the peaks with the new UB. " << std::endl;
    g_log.notice() << "Now, " << num_indexed << " are indexed with average error " << average_error << std::endl;
202

203
    // Save output properties
204

205
206
    this->setProperty("NumIndexed", num_indexed);
    this->setProperty("AverageError", average_error);
207
208
209
210
211
212
    }
  }


} // namespace Mantid
} // namespace Crystal