Skip to content

Commit

Permalink
Refs #12585 add integration
Browse files Browse the repository at this point in the history
  • Loading branch information
@VickieLynch
VickieLynch committed Jun 28, 2016
1 parent f49ff39 commit 4bda51f
Show file tree
Hide file tree
Showing 2 changed files with 111 additions and 39 deletions.
8 changes: 4 additions & 4 deletions Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegratePeaksMDHisto.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ class DLLExport IntegratePeaksMDHisto : public API::Algorithm {
public:
IntegratePeaksMDHisto();
/// Algorithm's name for identification
const std::string name() const override { return "IntegratePeaksMD"; };
const std::string name() const override { return "IntegratePeaksMDHisto"; };
/// Summary of algorithms purpose
const std::string summary() const override {
return "Integrate single-crystal peaks in reciprocal space, for "
Expand All @@ -42,9 +42,9 @@ class DLLExport IntegratePeaksMDHisto : public API::Algorithm {

DataObjects::MDHistoWorkspace_sptr normalize(
int h, int k, int l, double box, int gridPts,
API::MatrixWorkspace_const_sptr flux, API::MatrixWorkspace_const_sptr sa,
API::IMDEventWorkspace_const_sptr ws);

API::MatrixWorkspace_sptr flux, API::MatrixWorkspace_sptr sa,
API::IMDEventWorkspace_sptr ws);
void integratePeak(DataObjects::MDHistoWorkspace_sptr out, double& intensity, double& errorSquared, int gridPts);
};

} // namespace Mantid
Expand Down
142 changes: 107 additions & 35 deletions Framework/MDAlgorithms/src/IntegratePeaksMDHisto.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,14 @@
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VectorHelper.h"

#include "MantidGeometry/Crystal/IPeak.h"

namespace Mantid {
namespace MDAlgorithms {

using Mantid::Kernel::Direction;
using Mantid::API::WorkspaceProperty;
//using Mantid::API::WorkspaceProperty;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Kernel;
Expand Down Expand Up @@ -58,12 +57,6 @@ void IntegratePeaksMDHisto::init() {
"An input workspace containing momentum integrated vanadium "
"(a measure of the solid angle).");

declareProperty(make_unique<WorkspaceProperty<Workspace>>(
"OutputWorkspace", "", Direction::Output),
"A name for the output data MDHistoWorkspace.");
declareProperty(make_unique<WorkspaceProperty<Workspace>>(
"OutputNormalizationWorkspace", "", Direction::Output),
"A name for the output normalization MDHistoWorkspace.");
declareProperty(make_unique<WorkspaceProperty<PeaksWorkspace>>(
"PeaksWorkspace", "", Direction::Input),
"A PeaksWorkspace containing the peaks to integrate.");
Expand All @@ -81,77 +74,156 @@ void IntegratePeaksMDHisto::init() {
*/
void IntegratePeaksMDHisto::exec() {
/// Peak workspace to integrate
Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS =
PeaksWorkspace_sptr inPeakWS =
getProperty("PeaksWorkspace");

/// Output peaks workspace, create if needed
Mantid::DataObjects::PeaksWorkspace_sptr peakWS =
PeaksWorkspace_sptr peakWS =
getProperty("OutputWorkspace");
if (peakWS != inPeakWS)
peakWS = inPeakWS->clone();

API::MatrixWorkspace_const_sptr flux = getProperty("FluxWorkspace");
API::MatrixWorkspace_const_sptr sa =
MatrixWorkspace_sptr flux = getProperty("FluxWorkspace");
MatrixWorkspace_sptr sa =
getProperty("SolidAngleWorkspace");

API::IMDEventWorkspace_const_sptr m_inputWS = getProperty("InputWorkspace");
IMDEventWorkspace_sptr m_inputWS = getProperty("InputWorkspace");

double box = 0.5;
int gridPts = 201;

int npeaks = peakWS->getNumberPeaks();

auto prog = make_unique<API::Progress>(this, 0.3, 1.0, npeaks);
auto prog = make_unique<Progress>(this, 0.3, 1.0, npeaks);
PARALLEL_FOR1(peakWS)
for (int i = 0; i < npeaks; i++) {
PARALLEL_START_INTERUPT_REGION

IPeak &p = peakWS->getPeak(i);
// round to integer
int h = int(p.getH() + 0.5);
int k = int(p.getK() + 0.5);
int l = int(p.getL() + 0.5);
MDHistoWorkspace_sptr normBox = normalize(h, k, l, box, gridPts, flux, sa, m_inputWS);

//PARALLEL_CRITICAL(updateMD) {
//signal += m_normWS->getSignalAt(linIndex);
//}
int h = static_cast<int>(p.getH() + 0.5);
int k = static_cast<int>(p.getK() + 0.5);
int l = static_cast<int>(p.getL() + 0.5);
MDHistoWorkspace_sptr normBox = normalize(
h, k, l, box, gridPts, flux, sa, m_inputWS);
double intensity = 0.0;
double errorSquared = 0.0;
integratePeak(normBox, intensity, errorSquared, gridPts);
p.setIntensity(intensity);
p.setSigmaIntensity(sqrt(errorSquared));
prog->report();

PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}

MDHistoWorkspace_sptr IntegratePeaksMDHisto::normalize(
int h, int k, int l, double box, int gridPts,
API::MatrixWorkspace_const_sptr flux, API::MatrixWorkspace_const_sptr sa, API::IMDEventWorkspace_const_sptr ws) {

API::IAlgorithm_sptr normAlg = createChildAlgorithm("MDNormSCD");
MatrixWorkspace_sptr flux, MatrixWorkspace_sptr sa,
IMDEventWorkspace_sptr ws) {
IAlgorithm_sptr normAlg = createChildAlgorithm("MDNormSCD");
normAlg->setProperty("InputWorkspace", ws);
normAlg->setProperty("AlignedDim0",
"[H,0,0],"+boost::lexical_cast<std::string>(h-box)+","+boost::lexical_cast<std::string>(h+box)+","+boost::lexical_cast<std::string>(gridPts));
"[H,0,0],"+boost::lexical_cast<std::string>(h-box)+","+
boost::lexical_cast<std::string>(h+box)+","+
boost::lexical_cast<std::string>(gridPts));
normAlg->setProperty("AlignedDim1",
"[0,K,0],"+boost::lexical_cast<std::string>(k-box)+","+boost::lexical_cast<std::string>(k+box)+","+boost::lexical_cast<std::string>(gridPts));
"[0,K,0],"+boost::lexical_cast<std::string>(k-box)+","+
boost::lexical_cast<std::string>(k+box)+","+
boost::lexical_cast<std::string>(gridPts));
normAlg->setProperty("AlignedDim2",
"[0,0,L],"+boost::lexical_cast<std::string>(l-box)+","+boost::lexical_cast<std::string>(l+box)+","+boost::lexical_cast<std::string>(gridPts));
"[0,0,L],"+boost::lexical_cast<std::string>(l-box)+","+
boost::lexical_cast<std::string>(l+box)+","+
boost::lexical_cast<std::string>(gridPts));
normAlg->setProperty("FluxWorkspace", flux);
normAlg->setProperty("SolidAngleWorkspace", sa);
normAlg->setProperty("OutputWorkspace", "mdout");
normAlg->setProperty("OutputNormalizationWorkspace", "mdnorm");
normAlg->executeAsChildAlg();

MDHistoWorkspace_sptr mdout = normAlg->getProperty("OutputWorkspace");
MDHistoWorkspace_sptr mdnorm = normAlg->getProperty("OutputtNormalizationWorkspace");
Workspace_sptr mdout = normAlg->getProperty("OutputWorkspace");
Workspace_sptr mdnorm = normAlg->getProperty("OutputNormalizationWorkspace");

API::IAlgorithm_sptr alg = createChildAlgorithm("DivideMD");
IAlgorithm_sptr alg = createChildAlgorithm("DivideMD");
alg->setProperty("LHSWorkspace", mdout);
alg->setProperty("RHSWorkspace", mdnorm);
alg->setPropertyValue("OutputWorkspace", "out");
alg->execute();
MDHistoWorkspace_sptr out = alg->getProperty("OutputWorkspace");
return out;
Workspace_sptr out = alg->getProperty("OutputWorkspace");
return boost::dynamic_pointer_cast<MDHistoWorkspace>(out);
}

void IntegratePeaksMDHisto::integratePeak(MDHistoWorkspace_sptr out, double& intensity, double& errorSquared, int gridPts) {
double *F = out->getSignalArray();
int noPoints = 10;
double Fmax = 0;
double Fmin = 1e300;
for (int i = 0; i < gridPts*gridPts*gridPts; i++) {
if (!boost::math::isnan(F[i]) && !boost::math::isinf(F[i]) && F[i] != 0.0) {
if (F[i] < Fmin) Fmin = F[i];
if (F[i] > Fmax) Fmax = F[i];
}
}

double *SqError = out->getErrorSquaredArray();

double minIntensity = Fmin + 0.01 * (Fmax - Fmin);
int measuredPoints = 0;
int peakPoints = 0;
double peakSum = 0.0;
double measuredSum = 0.0;
double errSqSum = 0.0;
double measuredErrSqSum = 0.0;
//double Peak[201*201*201] = {0};
//double PeakErr2[201*201*201] = {0};
for (int Hindex = 0; Hindex < gridPts; Hindex++) {
for (int Kindex = 0; Kindex < gridPts; Kindex++) {
for (int Lindex = 0; Lindex < gridPts; Lindex++) {
int iHKL = Hindex + gridPts * (Kindex +gridPts * Lindex);
if (!boost::math::isnan(F[iHKL]) && !boost::math::isinf(F[iHKL])) {
measuredPoints = measuredPoints + 1;
measuredSum = measuredSum + F[iHKL];
measuredErrSqSum = measuredErrSqSum + F[iHKL];
if (F[iHKL] > minIntensity) {
int neighborPoints = 0;
for (int Hj = -2; Hj < 3; Hj++) {
for (int Kj = -2; Kj < 3; Kj++) {
for (int Lj = -2; Lj < 3; Lj++) {
int jHKL = Hindex + Hj + gridPts * (Kindex + Kj +gridPts * (Lindex + Lj));
if (Lindex+Lj >= 0 && Lindex+Lj < gridPts && Kindex+Kj >= 0 && Kindex+Kj < gridPts && Hindex+Hj >= 0 && Hindex+Hj < gridPts && F[jHKL] > minIntensity) {
neighborPoints = neighborPoints + 1;
}
}
}
}
if (neighborPoints >= noPoints) {
//Peak[iHKL] = F[iHKL];
//PeakErr2[iHKL] = SqError[iHKL];
peakPoints = peakPoints + 1;
peakSum = peakSum + F[iHKL];
errSqSum = errSqSum + SqError[iHKL];
}
}
}
else{
double minR = sqrt( std::pow(float(Hindex)/float(gridPts) - 0.5, 2) + std::pow(float(Kindex)/float(gridPts) - 0.5, 2) + std::pow(float(Lindex)/float(gridPts) - 0.5, 2));
if (minR < 0.1) {
intensity = 0.0;
errorSquared = 0.0;
return;
}
}
}
}
}
double ratio = float(peakPoints)/float(measuredPoints - peakPoints);
//std::cout peakSum, errSqSum, ratio, measuredSum, measuredErrSqSum
intensity = peakSum - ratio * (measuredSum - peakSum);
errorSquared = errSqSum + ratio * (measuredErrSqSum - errSqSum);
return;
}





Expand Down

0 comments on commit 4bda51f

Please sign in to comment.