Skip to content

Commit

Permalink
Refs #6338 PredictPeaks needs all goniometers for completeness
Browse files Browse the repository at this point in the history
  • Loading branch information
@VickieLynch
VickieLynch committed Dec 18, 2014
1 parent 8b980b7 commit 1d85aeb
Show file tree
Hide file tree
Showing 3 changed files with 160 additions and 112 deletions.
242 changes: 134 additions & 108 deletions Code/Mantid/Framework/Crystal/src/PredictPeaks.cpp
View file
Expand Up @@ -169,15 +169,51 @@ namespace Crystal
MatrixWorkspace_sptr matrixWS = boost::dynamic_pointer_cast<MatrixWorkspace>(inBareWS);
PeaksWorkspace_sptr peaksWS = boost::dynamic_pointer_cast<PeaksWorkspace>(inBareWS);
IMDEventWorkspace_sptr mdWS = boost::dynamic_pointer_cast<IMDEventWorkspace>(inBareWS);
std::vector<Matrix<double> > gonioVec;
gonio = Matrix<double>(3,3, true);
Mantid::Kernel::DblMatrix gonioLast = Matrix<double>(3,3, true);
if (matrixWS)
{
inWS = matrixWS;
// Retrieve the goniometer rotation matrix
try
{
gonio = inWS->mutableRun().getGoniometerMatrix();
gonioVec.push_back(gonio);
}
catch (std::runtime_error & e)
{
g_log.error() << "Error getting the goniometer rotation matrix from the InputWorkspace." << std::endl
<< e.what() << std::endl;
g_log.warning() << "Using identity goniometer rotation matrix instead." << std::endl;
}
}
else if (peaksWS)
{
inWS = peaksWS;
for (int i=0; i < static_cast<int>(peaksWS->getNumberPeaks()); ++i)
{

IPeak & p = peaksWS->getPeak(i);
gonio = p.getGoniometerMatrix();
if (!(gonio == gonioLast))
{
gonioLast = gonio;
gonioVec.push_back(gonioLast);
}
} // for each hkl in the workspace
}
else if (mdWS)
{
if (mdWS->getNumExperimentInfo() <= 0)
throw std::invalid_argument("Specified a MDEventWorkspace as InputWorkspace but it does not have any ExperimentInfo associated. Please choose a workspace with a full instrument and sample.");
inWS = mdWS->getExperimentInfo(0);
// Retrieve the goniometer rotation matrix
for (uint16_t i=0; i < mdWS->getNumExperimentInfo(); i++)
{
gonio =mdWS->getExperimentInfo(i)->mutableRun().getGoniometerMatrix();
gonioVec.push_back(gonio);
}
}
// Find the run number
runNumber = inWS->getRunNumber();
Expand Down Expand Up @@ -220,129 +256,119 @@ namespace Crystal
// Retrieve the OrientedLattice (UnitCell) from the workspace
crystal = inWS->sample().getOrientedLattice();

// Counter of possible peaks
numInRange = 0;

// Get the UB matrix from it
Matrix<double> ub(3,3, true);
ub = crystal.getUB();

// Retrieve the goniometer rotation matrix
gonio = Matrix<double>(3,3, true);
try
{
gonio = inWS->mutableRun().getGoniometerMatrix();
}
catch (std::runtime_error & e)
for (size_t iVec=0; iVec<gonioVec.size(); ++iVec)
{
g_log.error() << "Error getting the goniometer rotation matrix from the InputWorkspace." << std::endl
<< e.what() << std::endl;
g_log.warning() << "Using identity goniometer rotation matrix instead." << std::endl;
}
gonio = gonioVec[iVec];
// Final transformation matrix (HKL to Q in lab frame)
mat = gonio * ub;

// Final transformation matrix (HKL to Q in lab frame)
mat = gonio * ub;
// Sample position
V3D samplePos = inst->getSample()->getPos();

// Sample position
V3D samplePos = inst->getSample()->getPos();
// L1 path and direction
V3D beamDir = inst->getSource()->getPos() - samplePos;
//double L1 = beamDir.normalize(); // Normalize to unity

// L1 path and direction
V3D beamDir = inst->getSource()->getPos() - samplePos;
//double L1 = beamDir.normalize(); // Normalize to unity
if ((fabs(beamDir.X()) > 1e-2) || (fabs(beamDir.Y()) > 1e-2)) // || (beamDir.Z() < 0))
throw std::invalid_argument("Instrument must have a beam direction that is only in the +Z direction for this algorithm to be valid..");

if ((fabs(beamDir.X()) > 1e-2) || (fabs(beamDir.Y()) > 1e-2)) // || (beamDir.Z() < 0))
throw std::invalid_argument("Instrument must have a beam direction that is only in the +Z direction for this algorithm to be valid..");

// Counter of possible peaks
numInRange = 0;

if (HKLPeaksWorkspace)
{
// --------------Use the HKL from a list in a PeaksWorkspace --------------------------
// Disable some of the other filters
minD = 0.0;
maxD = 1e10;
wlMin = 0.0;
wlMax = 1e10;

// PRAGMA_OMP(parallel for schedule(dynamic, 1) )
for (int i=0; i < static_cast<int>(HKLPeaksWorkspace->getNumberPeaks()); ++i)
{
PARALLEL_START_INTERUPT_REGION

IPeak & p = HKLPeaksWorkspace->getPeak(i);
// Get HKL from that peak
V3D hkl = p.getHKL();
// Use the rounded HKL value on option
if (RoundHKL)
hkl.round();
// Predict the HKL of that peak
doHKL(hkl[0], hkl[1], hkl[2], false);

PARALLEL_END_INTERUPT_REGION
} // for each hkl in the workspace
PARALLEL_CHECK_INTERUPT_REGION
}
else
{
// ---------------- Determine which HKL to look for -------------------------------------
// Inverse d-spacing that is the limit to look for.
double Qmax = 2.*M_PI/minD;
V3D hklMin(0,0,0);
V3D hklMax(0,0,0);
for (double qx=-1; qx < 2; qx += 2)
for (double qy=-1; qy < 2; qy += 2)
for (double qz=-1; qz < 2; qz += 2)
if (HKLPeaksWorkspace)
{
// --------------Use the HKL from a list in a PeaksWorkspace --------------------------
// Disable some of the other filters
minD = 0.0;
maxD = 1e10;
wlMin = 0.0;
wlMax = 1e10;

// PRAGMA_OMP(parallel for schedule(dynamic, 1) )
for (int i=0; i < static_cast<int>(HKLPeaksWorkspace->getNumberPeaks()); ++i)
{
// Build a q-vector for this corner of a cube
V3D Q(qx,qy,qz);
Q *= Qmax;
V3D hkl = crystal.hklFromQ(Q);
// Find the limits of each hkl
for (size_t i=0; i<3; i++)
{
if (hkl[i] < hklMin[i]) hklMin[i] = hkl[i];
if (hkl[i] > hklMax[i]) hklMax[i] = hkl[i];
}
}
// Round to nearest int
hklMin.round();
hklMax.round();

// How many HKLs is that total?
V3D hklDiff = hklMax-hklMin + V3D(1,1,1);
size_t numHKLs = size_t( hklDiff[0] * hklDiff[1] * hklDiff[2]);

g_log.information() << "HKL range for d_min of " << minD << "to d_max of " << maxD << " is from " << hklMin << " to " << hklMax << ", a total of " << numHKLs << " possible HKL's\n";

if (numHKLs > 10000000000)
throw std::invalid_argument("More than 10 billion HKLs to search. Is your d_min value too small?");

Progress prog(this, 0.0, 1.0, numHKLs);
prog.setNotifyStep(0.01);

// PRAGMA_OMP(parallel for schedule(dynamic, 1) )
for (int h=(int)hklMin[0]; h <= (int)hklMax[0]; h++)
{
PARALLEL_START_INTERUPT_REGION
for (int k=(int)hklMin[1]; k <= (int)hklMax[1]; k++)
PARALLEL_START_INTERUPT_REGION

IPeak & p = HKLPeaksWorkspace->getPeak(i);
// Get HKL from that peak
V3D hkl = p.getHKL();
// Use the rounded HKL value on option
if (RoundHKL)
hkl.round();
// Predict the HKL of that peak
doHKL(hkl[0], hkl[1], hkl[2], false);

PARALLEL_END_INTERUPT_REGION
} // for each hkl in the workspace
PARALLEL_CHECK_INTERUPT_REGION
}
else
{
for (int l=(int)hklMin[2]; l <= (int)hklMax[2]; l++)
// ---------------- Determine which HKL to look for -------------------------------------
// Inverse d-spacing that is the limit to look for.
double Qmax = 2.*M_PI/minD;
V3D hklMin(0,0,0);
V3D hklMax(0,0,0);
for (double qx=-1; qx < 2; qx += 2)
for (double qy=-1; qy < 2; qy += 2)
for (double qz=-1; qz < 2; qz += 2)
{
// Build a q-vector for this corner of a cube
V3D Q(qx,qy,qz);
Q *= Qmax;
V3D hkl = crystal.hklFromQ(Q);
// Find the limits of each hkl
for (size_t i=0; i<3; i++)
{
if (hkl[i] < hklMin[i]) hklMin[i] = hkl[i];
if (hkl[i] > hklMax[i]) hklMax[i] = hkl[i];
}
}
// Round to nearest int
hklMin.round();
hklMax.round();

// How many HKLs is that total?
V3D hklDiff = hklMax-hklMin + V3D(1,1,1);
size_t numHKLs = size_t( hklDiff[0] * hklDiff[1] * hklDiff[2]);

g_log.information() << "HKL range for d_min of " << minD << "to d_max of " << maxD << " is from " << hklMin << " to " << hklMax << ", a total of " << numHKLs << " possible HKL's\n";

if (numHKLs > 10000000000)
throw std::invalid_argument("More than 10 billion HKLs to search. Is your d_min value too small?");

Progress prog(this, 0.0, 1.0, numHKLs);
prog.setNotifyStep(0.01);

// PRAGMA_OMP(parallel for schedule(dynamic, 1) )
for (int h=(int)hklMin[0]; h <= (int)hklMax[0]; h++)
{
if (refCond->isAllowed(h,k,l) && (abs(h) + abs(k) + abs(l) != 0))
PARALLEL_START_INTERUPT_REGION
for (int k=(int)hklMin[1]; k <= (int)hklMax[1]; k++)
{
doHKL(double(h), double(k), double(l), true);
} // refl is allowed and not 0,0,0
prog.report();
} // for each l
} // for each k
PARALLEL_END_INTERUPT_REGION
} // for each h
PARALLEL_CHECK_INTERUPT_REGION

} // Find the HKL automatically

g_log.notice() << "Out of " << numInRange << " allowed peaks within parameters, " << pw->getNumberPeaks() << " were found to hit a detector.\n";
for (int l=(int)hklMin[2]; l <= (int)hklMax[2]; l++)
{
if (refCond->isAllowed(h,k,l) && (abs(h) + abs(k) + abs(l) != 0))
{
doHKL(double(h), double(k), double(l), true);
} // refl is allowed and not 0,0,0
prog.report();
} // for each l
} // for each k
PARALLEL_END_INTERUPT_REGION
} // for each h
PARALLEL_CHECK_INTERUPT_REGION

} // Find the HKL automatically

g_log.notice() << "Out of " << numInRange << " allowed peaks within parameters, " << pw->getNumberPeaks() << " were found to hit a detector.\n";
}
}



} // namespace Mantid
} // namespace Crystal
21 changes: 18 additions & 3 deletions Code/Mantid/Framework/Crystal/src/SortHKL.cpp
View file
Expand Up @@ -91,6 +91,7 @@ namespace Crystal
tablews->addColumn("double", "Rmerge");
tablews->addColumn("double", "Rpim");
tablews->addColumn("double", "Data Completeness");

// append to the table workspace
API::TableRow newrow = tablews->appendRow();
newrow << "Overall";
Expand Down Expand Up @@ -134,13 +135,26 @@ namespace Crystal
}
}
std::vector< std::pair<std::string, bool> > criteria;
// Sort by bank and then HKL
criteria.push_back( std::pair<std::string, bool>("Wavelength", true) );
// Sort by wavelength
criteria.push_back( std::pair<std::string, bool>("wavelength", true) );
peaksW->sort(criteria);
int unique = NumberPeaks - equivalent;
// table workspace output
newrow << unique << peaks[0].getWavelength() << peaks[NumberPeaks-1].getWavelength();

API::IAlgorithm_sptr predictAlg = createChildAlgorithm("PredictPeaks");
predictAlg->setProperty("InputWorkspace", InPeaksW);
predictAlg->setPropertyValue("OutputWorkspace", "predictedPeaks");
predictAlg->setProperty("WavelengthMin",peaks[0].getWavelength());
predictAlg->setProperty("WavelengthMax",peaks[NumberPeaks-1].getWavelength());
// Sort by dspacing
criteria.push_back( std::pair<std::string, bool>("dspacing", true) );
peaksW->sort(criteria);
predictAlg->setProperty("MinDSpacing",peaks[0].getDSpacing());
predictAlg->executeAsChildAlg();
PeaksWorkspace_sptr predictedWksp = predictAlg->getProperty("OutputWorkspace");
int predictedPeaks = predictedWksp->getNumberPeaks();

criteria.clear();
// Sort by bank and then HKL
//criteria.push_back( std::pair<std::string, bool>("BankName", true) );
Expand Down Expand Up @@ -240,7 +254,8 @@ namespace Crystal
Statistics statsMult = getStatistics(multiplicity);
multiplicity.clear();
// statistics to output table workspace
newrow << statsMult.mean << statsIsigI.mean << 100.0 * rSum / f2Sum << 100.0 * rpSum / f2Sum << 0.0;
newrow << statsMult.mean << statsIsigI.mean << 100.0 * rSum / f2Sum << 100.0 * rpSum / f2Sum
<< double(NumberPeaks)/double(predictedPeaks);
data.clear();
sig2.clear();
//Reset hkl of equivalent peaks to original value
Expand Down
9 changes: 8 additions & 1 deletion Code/Mantid/Framework/Crystal/test/SortHKLTest.h
View file
Expand Up @@ -8,6 +8,7 @@
#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>
Expand Down Expand Up @@ -38,7 +39,13 @@ class SortHKLTest : public CxxTest::TestSuite
Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentRectangular(4, 10, 1.0);
PeaksWorkspace_sptr ws(new PeaksWorkspace());
ws->setInstrument(inst);
//ws->setName("TOPAZ_peaks");

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

double smu = 0.357;
double amu = 0.011;
NeutronAtom neutron(static_cast<uint16_t>(EMPTY_DBL()), static_cast<uint16_t>(0),
Expand Down

0 comments on commit 1d85aeb

Please sign in to comment.