Skip to content

Commit

Permalink
Refs #6338 options for resolution shells and banks plus tests
Browse files Browse the repository at this point in the history
  • Loading branch information
VickieLynch committed Jan 16, 2015
1 parent 84fda71 commit 516d3e0
Show file tree
Hide file tree
Showing 5 changed files with 226 additions and 7 deletions.
1 change: 1 addition & 0 deletions Code/Mantid/Framework/Crystal/CMakeLists.txt
Expand Up @@ -202,6 +202,7 @@ set ( TEST_FILES
ShowPossibleCellsTest.h
SortHKLTest.h
SortPeaksWorkspaceTest.h
StatisticsOfPeaksWorkspaceTest.h
TransformHKLTest.h
)

Expand Down
5 changes: 5 additions & 0 deletions Code/Mantid/Framework/Crystal/src/SortHKL.cpp
Expand Up @@ -125,6 +125,11 @@ void SortHKL::exec() {
peaksW->removePeak(i);
}
NumberPeaks = peaksW->getNumberPeaks();
if (NumberPeaks == 0)
{
g_log.error() << "Number of peaks should not be 0 for SortHKL.\n";
return;
}
int equivalent = 0;
for (int i = 0; i < NumberPeaks; i++) {
V3D hkl1 = peaks[i].getHKL();
Expand Down
55 changes: 48 additions & 7 deletions Code/Mantid/Framework/Crystal/src/StatisticsOfPeaksWorkspace.cpp
Expand Up @@ -60,8 +60,16 @@ void StatisticsOfPeaksWorkspace::init() {
Direction::Output),
"Output PeaksWorkspace");
declareProperty(new WorkspaceProperty<ITableWorkspace>(
"StatisticsTable", "StaticticsTable", Direction::Output),
"StatisticsTable", "StatisticsTable", Direction::Output),
"An output table workspace for the statistics of the peaks.");
std::vector<std::string> sortTypes;
sortTypes.push_back("ResolutionShell");
sortTypes.push_back("Bank");
sortTypes.push_back("RunNumber");
sortTypes.push_back("Overall");
declareProperty("SortBy", sortTypes[0],
boost::make_shared<StringListValidator>(sortTypes),
"Sort the peaks by bank, run number(default) or only overall statistics.");
}

//----------------------------------------------------------------------------------------------
Expand All @@ -70,13 +78,17 @@ void StatisticsOfPeaksWorkspace::init() {
void StatisticsOfPeaksWorkspace::exec() {

ws = getProperty("InputWorkspace");
std::string sortType = getProperty("SortBy");
IPeaksWorkspace_sptr tempWS = WorkspaceFactory::Instance().createPeaks();
// Copy over ExperimentInfo from input workspace
tempWS->copyExperimentInfoFrom(ws.get());

// We must sort the peaks
std::vector<std::pair<std::string, bool>> criteria;
criteria.push_back(std::pair<std::string, bool>("RunNumber", true));
if (sortType.compare(0, 2, "Re") == 0)
criteria.push_back(std::pair<std::string, bool>("wavelength", false));
else if (sortType.compare(0, 2, "Ru") == 0)
criteria.push_back(std::pair<std::string, bool>("RunNumber", true));
criteria.push_back(std::pair<std::string, bool>("BankName", true));
criteria.push_back(std::pair<std::string, bool>("h", true));
criteria.push_back(std::pair<std::string, bool>("k", true));
Expand All @@ -88,24 +100,53 @@ void StatisticsOfPeaksWorkspace::exec() {
// =========================================
std::vector<Peak> peaks = ws->getPeaks();
doSortHKL(ws, "Overall");
if (sortType.compare(0, 2, "Ov") == 0) return;
// run number
int oldSequence = peaks[0].getRunNumber();
std::string oldSequence;
if (sortType.compare(0, 2, "Re") == 0)
{
double wavelength = peaks[0].getWavelength();
if (wavelength > 3.0) oldSequence = "first";
else if (wavelength > 2.5) oldSequence = "second";
else if (wavelength > 2.0) oldSequence = "third";
else if (wavelength > 1.5) oldSequence = "fourth";
else if (wavelength > 1.0) oldSequence = "fifth";
else if (wavelength > 0.5) oldSequence = "sixth";
else oldSequence = "seventh";
}
else if (sortType.compare(0, 2, "Ru") == 0)
oldSequence = boost::lexical_cast<std::string>(peaks[0].getRunNumber());
else oldSequence= peaks[0].getBankName();
// Go through each peak at this run / bank
for (int wi = 0; wi < ws->getNumberPeaks(); wi++) {

Peak &p = peaks[wi];
int sequence = p.getRunNumber();
std::string sequence;
if (sortType.compare(0, 2, "Re") == 0)
{
double wavelength = p.getWavelength();
if (wavelength > 3.0) sequence = "first";
else if (wavelength > 2.5) sequence = "second";
else if (wavelength > 2.0) sequence = "third";
else if (wavelength > 1.5) sequence = "fourth";
else if (wavelength > 1.0) sequence = "fifth";
else if (wavelength > 0.5) sequence = "sixth";
else sequence = "seventh";
}
else if (sortType.compare(0, 2, "Ru") == 0)
sequence = boost::lexical_cast<std::string>(p.getRunNumber());
else sequence= p.getBankName();

if (sequence != oldSequence && tempWS->getNumberPeaks()>0) {
doSortHKL(tempWS, boost::lexical_cast<std::string>(oldSequence));
if (sequence.compare(oldSequence) !=0 && tempWS->getNumberPeaks()>0) {
doSortHKL(tempWS, oldSequence);
tempWS = WorkspaceFactory::Instance().createPeaks();
// Copy over ExperimentInfo from input workspace
tempWS->copyExperimentInfoFrom(ws.get());
oldSequence = sequence;
}
tempWS->addPeak(p);
}
doSortHKL(tempWS, boost::lexical_cast<std::string>(oldSequence));
doSortHKL(tempWS, oldSequence);
}
//----------------------------------------------------------------------------------------------
/** Perform SortHKL on the output workspaces
Expand Down
119 changes: 119 additions & 0 deletions Code/Mantid/Framework/Crystal/test/StatisticsOfPeaksWorkspaceTest.h
@@ -0,0 +1,119 @@
#ifndef MANTID_CRYSTAL_STATISTICSOFPEAKSWORKSPACETEST_H_
#define MANTID_CRYSTAL_STATISTICSOFPEAKSWORKSPACETEST_H_

#include "MantidCrystal/StatisticsOfPeaksWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/IDTypes.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Timer.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include <cxxtest/TestSuite.h>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <Poco/File.h>

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

class StatisticsOfPeaksWorkspaceTest : public CxxTest::TestSuite
{
public:

void test_Init()
{
StatisticsOfPeaksWorkspace alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
}

void do_test(int numRuns, size_t numBanks, size_t numPeaksPerBank)
{
Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentRectangular(4, 10, 1.0);
PeaksWorkspace_sptr ws(new PeaksWorkspace());
ws->setInstrument(inst);

auto lattice = new Mantid::Geometry::OrientedLattice;
Mantid::Kernel::DblMatrix UB(3, 3, true);
UB.identityMatrix();
lattice->setUB(UB);
ws->mutableSample().setOrientedLattice(lattice);

for (int run=1000; run<numRuns+1000; run++)
for (size_t b=1; b<=numBanks; b++)
for (size_t i=0; i<numPeaksPerBank; i++)
{
V3D hkl(static_cast<double>(i), static_cast<double>(i), static_cast<double>(i));
DblMatrix gon(3,3, true);
Peak p(inst, static_cast<detid_t>(b*100 + i+1+i*10), static_cast<double>(i)*1.0+0.5, hkl, gon);
p.setRunNumber(run);
p.setBankName("bank1");
p.setIntensity( static_cast<double>(i)+0.1);
p.setSigmaIntensity( sqrt(static_cast<double>(i)+0.1));
p.setBinCount( static_cast<double>(i) );
ws->addPeak(p);
}
AnalysisDataService::Instance().addOrReplace("TOPAZ_peaks", ws);

StatisticsOfPeaksWorkspace alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setProperty("InputWorkspace", "TOPAZ_peaks") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("SortBy", "Overall") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("StatisticsTable", "stat") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("OutputWorkspace", "TOPAZ_peaks") );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );

ITableWorkspace_sptr tableOut;
TS_ASSERT_THROWS_NOTHING( tableOut = boost::dynamic_pointer_cast<ITableWorkspace>(
AnalysisDataService::Instance().retrieve("stat") ) );
TS_ASSERT(tableOut);
if (!tableOut) return;
TS_ASSERT_EQUALS( tableOut->rowCount(),1);
TS_ASSERT_EQUALS( tableOut->String(0,0), "Overall");
TS_ASSERT_EQUALS( tableOut->Int(0,1), 3);
TS_ASSERT_DELTA( tableOut->Double(0,2), 1.5,.1);
TS_ASSERT_DELTA( tableOut->Double(0,3), 3.5,.1);
TS_ASSERT_DELTA( tableOut->Double(0,4), 6.,.1);
TS_ASSERT_DELTA( tableOut->Double(0,5), 1.4195,.1);
TS_ASSERT_DELTA( tableOut->Double(0,6), 0.,.1);
TS_ASSERT_DELTA( tableOut->Double(0,7), 0.,.1);

PeaksWorkspace_sptr wsout;
TS_ASSERT_THROWS_NOTHING( wsout = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve("TOPAZ_peaks") ) );
TS_ASSERT(wsout);
if (!wsout) return;
TS_ASSERT_EQUALS( wsout->getNumberPeaks(), 24);

Peak p = wsout->getPeaks()[0];
TS_ASSERT_EQUALS(p.getH(),1 );
TS_ASSERT_EQUALS(p.getK(),1 );
TS_ASSERT_EQUALS(p.getL(),1 );
TS_ASSERT_DELTA(p.getIntensity(),1.1, 1e-4);
TS_ASSERT_DELTA(p.getSigmaIntensity(),1.0488, 1e-4);
TS_ASSERT_DELTA(p.getWavelength(),1.5, 1e-4 );
TS_ASSERT_EQUALS(p.getRunNumber(),1000. );
TS_ASSERT_DELTA(p.getDSpacing(),3.5933, 1e-4 );
}

/// Test with a few peaks
void test_exec()
{
do_test(2, 4,4);
}


};


#endif /* MANTID_CRYSTAL_STATISTICSOFPEAKSWORKSPACETEST_H_ */

@@ -0,0 +1,53 @@
.. algorithm::

.. summary::

.. alias::

.. properties::

Description
-----------

Statistics of the Peaks Workspaces are calculated for all peaks and by
default for resolution shell. There is a SortBy option to change this
to by orientation (RunNumber) or by Anger camera (bank) or only do all peaks.

Statistics include:
No of Unique Reflections
Resolution range
Multiplicity
Mean ((I)/sd(I))
Rmerge
Rplm
Data Completeness


Usage
-----

**Example - an example of running StatisticsOfPeaksWorkspace with PointGroup option.**

.. testcode:: ExStatisticsOfPeaksWorkspaceOption

#load a peaks workspace from file
peaks = LoadIsawPeaks(Filename=r'Peaks5637.integrate')
LoadIsawUB(peaks,"ls5637.mat")
peak = peaks.getPeak(0)
print "HKL of first peak in table %d" % peak.getH(),peak.getK(),peak.getL()

pw,stats = StatisticsOfPeaksWorkspace(InputWorkspace=peaks, PointGroup='-31m (Trigonal - Rhombohedral)', SortBy="Overall")
peak = pw.getPeak(0)
print "HKL of first peak in table %d" % peak.getH(),peak.getK(),peak.getL()
print stats.row(0)

Output:

.. testoutput:: ExStatisticsOfPeaksWorkspaceOption

HKL of first peak in table 1 4.0 -9.0
HKL of first peak in table -10 3.0 -40.0
{'Data Completeness': inf, 'Rmerge': 76.453114067956719, 'Multiplicity': 1.0074074074074073, 'Resolution Min': 0.29574100000000003, 'No. of Unique Reflections': 405, 'Mean ((I)/sd(I))': 27.50726166791943, 'Resolution Max': 3.1660760000000003, 'Resolution Shell': 'Overall', 'Rpim': 76.453114067956719}


.. categories::

0 comments on commit 516d3e0

Please sign in to comment.