Skip to content

Commit

Permalink
RE 5672, Added option to round or not round hkl to IndexPeaks()
Browse files Browse the repository at this point in the history 
Added parameter RoundHKLs to IndexPeaks() algorithm, and updated
IndexPeaksTest to verify that the peaks are rounded if RoundHKLs
is set true and are NOT rounded if RoundHKLs is set false.
refs #5672
  • Loading branch information
@DennisMikkelson
DennisMikkelson committed Aug 2, 2012
1 parent baffaee commit 3bdd99c
Show file tree
Hide file tree
Showing 2 changed files with 96 additions and 14 deletions.
35 changes: 30 additions & 5 deletions Code/Mantid/Framework/Crystal/src/IndexPeaks.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,17 @@
/*WIKI*
Given a PeaksWorkspace with a UB matrix stored with the sample, this algorithm will use UB inverse to index the peaks.
Any peak with any Miller index more than the specified tolerance away from an integer will have its (h,k,l) set to (0,0,0).
Given a PeaksWorkspace with a UB matrix stored with the sample, this algorithm will use UB inverse
to index the peaks. If there are peaks from multiple runs in the workspace, the stored UB will be
used to get an initial indexing for the peaks from each individual run. Subsequently, a temporary
UB will be optimzed for the peaks from each individual run, and used to index the peaks from that
run. In this way, a consistent indexing of the peaks from multiple runs will be obtained, which
indexes the largest number of peaks, although one UB may not produce exactly that indexing for
all peaks, within the specified tolerance.
Any peak with any Miller index more than the specified tolerance away from an integer will have its
(h,k,l) set to (0,0,0). The calculated Miller indices can either be rounded to the nearest integer
value, or can be left as decimal fractions.
*WIKI*/
Expand Down Expand Up @@ -72,6 +81,8 @@ namespace Crystal

this->declareProperty(new 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");
}

//--------------------------------------------------------------------------
Expand All @@ -98,6 +109,8 @@ namespace Crystal
"ERROR: The stored UB is not a valid orientation matrix");
}

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

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

Expand Down Expand Up @@ -147,7 +160,9 @@ namespace Crystal
tolerance,
miller_indices,
original_error );
IndexingUtils::RoundHKLs( miller_indices );

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

Expand All @@ -174,9 +189,11 @@ namespace Crystal
tolerance,
miller_indices,
average_error );
IndexingUtils::RoundHKLs( miller_indices );

if ( num_indexed < original_indexed ) // just use the original UB
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;
Expand All @@ -186,6 +203,14 @@ namespace Crystal
iteration++;
}

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

total_indexed += num_indexed;
total_error = average_error * num_indexed;

Expand Down
75 changes: 66 additions & 9 deletions Code/Mantid/Framework/Crystal/test/IndexPeaksTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ class IndexPeaksTest : public CxxTest::TestSuite
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", WSName) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Tolerance","0.1") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RoundHKLs",false) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );

Expand All @@ -89,30 +90,86 @@ class IndexPeaksTest : public CxxTest::TestSuite
// Check the output properties
int numIndexed = alg.getProperty("NumIndexed");
TS_ASSERT_EQUALS( numIndexed, 43);

double averageError = alg.getProperty("AverageError");
TS_ASSERT_DELTA( averageError, 0.0097, 1e-3);

// spot check a few peaks
V3D peak_0_hkl ( -4, -1, -6 ); // first peak
// spot check a few peaks for
// fractional Miller indices

V3D peak_0_hkl_d ( -4.03065, -0.9885090, -6.01095 ); // first peak
V3D peak_1_hkl_d ( -2.99276, 0.9955220, -4.00375 );
V3D peak_2_hkl_d ( -3.99311, 0.9856010, -5.00772 );
V3D peak_10_hkl_d( -3.01107, -0.0155531, -7.01377 );
V3D peak_42_hkl_d( -1.97065, 4.0283600, -6.97828 ); // last peak

V3D error = peak_0_hkl_d -peaks[ 0].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-4 );

error = peak_1_hkl_d -peaks[ 1].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-4 );

error = peak_2_hkl_d -peaks[ 2].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-4 );

error = peak_10_hkl_d -peaks[10].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-4 );

error = peak_42_hkl_d -peaks[42].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-4 );

// clear all the peak indexes then
// re-run the algorithm, rounding the HKLs
// this time, and again check a few peaks
for ( size_t i = 0; i < n_peaks; i++ )
{
peaks[i].setHKL( V3D(0,0,0) );
}

TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", WSName) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Tolerance","0.1") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RoundHKLs",true) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );

// Check that the peaks were all indexed
for ( size_t i = 0; i < n_peaks; i++ )
{
TS_ASSERT( IndexingUtils::ValidIndex( peaks[i].getHKL(), tolerance ) );
}

// Check the output properties
numIndexed = alg.getProperty("NumIndexed");
TS_ASSERT_EQUALS( numIndexed, 43);

averageError = alg.getProperty("AverageError");
TS_ASSERT_DELTA( averageError, 0.0097, 1e-3);

// spot check a few peaks for
// integer Miller indices
V3D peak_0_hkl ( -4, -1, -6 );
V3D peak_1_hkl ( -3, 1, -4 );
V3D peak_2_hkl ( -4, 1, -5 );
V3D peak_10_hkl( -3, 0, -7 );
V3D peak_42_hkl( -2, 4, -7 ); // last peak
V3D error = peak_0_hkl -peaks[ 0].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 0.1 );

error = peak_0_hkl -peaks[ 0].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-10 );

error = peak_1_hkl -peaks[ 1].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 0.1 );
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-10 );

error = peak_2_hkl -peaks[ 2].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 0.1 );
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-10 );

error = peak_10_hkl -peaks[10].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 0.1 );
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-10 );

error = peak_42_hkl -peaks[42].getHKL();
TS_ASSERT_DELTA( error.norm(), 0.0, 0.1 );
TS_ASSERT_DELTA( error.norm(), 0.0, 1e-10 );

// Remove workspace from the data service.
AnalysisDataService::Instance().remove(WSName);
}
Expand Down

0 comments on commit 3bdd99c

Please sign in to comment.