IndexPeaks.cpp 7.47 KB
Newer Older
1
2
3
4
5
#include "MantidCrystal/IndexPeaks.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
6
#include "MantidKernel/BoundedValidator.h"
7
#include <cstdio>
8
9
10
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

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;

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

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

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

45
    boost::shared_ptr<BoundedValidator<double> > mustBePositive(new BoundedValidator<double>());
46
47
    mustBePositive->setLower(0.0);

48
    this->declareProperty(new PropertyWithValue<double>( "Tolerance",0.15,
49
          mustBePositive,Direction::Input),"Indexing Tolerance (0.15)");
50
51
52

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

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

    this->declareProperty( "RoundHKLs", true, "Round H, K and L values to integers");
58
59
60
61
62
63
64
  }

  //--------------------------------------------------------------------------
  /** Execute the algorithm.
   */
  void IndexPeaks::exec()
  {
65
    PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
66
67
68
69
70
71
72
73
    if (!ws) 
    { 
      throw std::runtime_error("Could not read the peaks workspace");
    }

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

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

80
81
    bool round_hkls = this->getProperty("RoundHKLs");

82
83
84
    std::vector<Peak> &peaks = ws->getPeaks();
    size_t n_peaks = ws->getNumberPeaks();

85
86
    // get list of run numbers in this peaks workspace
    std::vector<int> run_numbers;
87
88
    for ( size_t i = 0; i < n_peaks; i++ )
    {
89
90
91
92
93
94
95
96
97
98
99
100
      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 );
101
102
    }

103
104
105
    // index the peaks for each run separately, using a UB matrix optimized for
    // that run
    int    total_indexed = 0;
106
    double average_error;
107
108
    double total_error = 0;
    double tolerance = this->getProperty("Tolerance");
109

110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    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( tempUB, q_vectors, 
                                                                tolerance, 
                                                                miller_indices,
                                                                original_error );
131
132
133

      IndexingUtils::RoundHKLs( miller_indices );  // HKLs must be rounded for 
                                                   // Optimize_UB to work
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
      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( tempUB, q_vectors, 
157
158
159
                                                             tolerance, 
                                                             miller_indices,
                                                             average_error );
160

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

        if ( num_indexed < original_indexed )       // just use the original UB
165
166
167
168
169
170
171
172
173
        {
          num_indexed = original_indexed;
          average_error = original_error;
          done = true;
        }

        iteration++;
      }

174
175
176
177
178
179
180
181
      if ( !round_hkls )      // If user wants fractional hkls, recalculate them
      {
        num_indexed = IndexingUtils::CalculateMillerIndices( tempUB, q_vectors,
                                                             tolerance,
                                                             miller_indices,
                                                             average_error );
      }

182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
      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 << std::endl;
        g_log.notice() << "Average error in h,k,l for indexed peaks =  " << average_error << std::endl;
      }

      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++;
        }
      }
    }
203

204
205
206
207
    if ( total_indexed > 0 )
      average_error = total_error / total_indexed;
    else
      average_error = 0;
208

209
210
211
212
    // 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 << std::endl;
    g_log.notice() << "Average error in h,k,l for indexed peaks =  " << average_error << std::endl;
213
214

    // Save output properties
215
    this->setProperty("NumIndexed", total_indexed);
216
    this->setProperty("AverageError", average_error);
Lynch, Vickie's avatar
Lynch, Vickie committed
217
218
    // Show the lattice parameters
    g_log.notice() << o_lattice << "\n";
219
220
221
222
223
  }


} // namespace Mantid
} // namespace Crystal