-
Dimitar Tasev authoredDimitar Tasev authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
CrossCorrelate.cpp 9.73 KiB
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/CrossCorrelate.h"
#include "MantidAPI/RawCountValidator.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/VectorHelper.h"
#include <boost/iterator/counting_iterator.hpp>
#include <numeric>
#include <sstream>
namespace Mantid {
namespace Algorithms {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(CrossCorrelate)
using namespace Kernel;
using namespace API;
/// Initialisation method.
void CrossCorrelate::init() {
auto wsValidator = boost::make_shared<CompositeValidator>();
wsValidator->add<API::WorkspaceUnitValidator>("dSpacing");
wsValidator->add<API::RawCountValidator>();
// Input and output workspaces
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"InputWorkspace", "", Direction::Input, wsValidator),
"A 2D workspace with X values of d-spacing");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"The name of the output workspace");
auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
mustBePositive->setLower(0);
// Reference spectra against which cross correlation is performed
declareProperty("ReferenceSpectra", 0, mustBePositive,
"The Workspace Index of the spectra to correlate all other "
"spectra against. ");
// Spectra in the range [min to max] will be cross correlated to reference.
declareProperty("WorkspaceIndexMin", 0, mustBePositive,
"The workspace index of the first member of the range of "
"spectra to cross-correlate against.");
declareProperty("WorkspaceIndexMax", 0, mustBePositive,
" The workspace index of the last member of the range of "
"spectra to cross-correlate against.");
// Only the data in the range X_min, X_max will be used
declareProperty("XMin", 0.0,
"The starting point of the region to be cross correlated.");
declareProperty("XMax", 0.0,
"The ending point of the region to be cross correlated.");
}
/** Executes the algorithm
*
* @throw runtime_error Thrown if algorithm cannot execute
*/
void CrossCorrelate::exec() {
MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
int reference = getProperty("ReferenceSpectra");
const size_t index_ref = static_cast<size_t>(reference);
// check that the data range specified makes sense
double xmin = getProperty("XMin"); double xmax = getProperty("XMax");
if (xmin >= xmax)
throw std::runtime_error("Must specify xmin < xmax");
// Now check if the range between x_min and x_max is valid
auto &referenceX = inputWS->x(index_ref);
auto minIt = std::find_if(referenceX.cbegin(), referenceX.cend(),
std::bind2nd(std::greater<double>(), xmin));
if (minIt == referenceX.cend())
throw std::runtime_error("No data above XMin");
auto maxIt = std::find_if(minIt, referenceX.cend(),
std::bind2nd(std::greater<double>(), xmax));
if (minIt == maxIt)
throw std::runtime_error("Range is not valid");
MantidVec::difference_type difminIt =
std::distance(referenceX.cbegin(), minIt);
MantidVec::difference_type difmaxIt =
std::distance(referenceX.cbegin(), maxIt);
// Now loop on the spectra in the range spectra_min and spectra_max and get
// valid spectra
int specmin = getProperty("WorkspaceIndexMin");
int specmax = getProperty("WorkspaceIndexMax");
if (specmin >= specmax)
throw std::runtime_error(
"Must specify WorkspaceIndexMin<WorkspaceIndexMax");
// Get the number of spectra in range specmin to specmax
int nspecs = 1 + specmax - specmin;
// Indexes of all spectra in range
std::vector<size_t> indexes(boost::make_counting_iterator(specmin),
boost::make_counting_iterator(specmax + 1));
std::ostringstream mess;
if (nspecs == 0) // Throw if no spectra in range
{
mess << "No Workspaces in range between" << specmin << " and " << specmax;
throw std::runtime_error(mess.str());
}
// Output message information
mess << "There are " << nspecs << " Workspaces in the range\n";
g_log.information(mess.str());
mess.str("");
// Take a copy of the reference spectrum
auto &referenceY = inputWS->y(index_ref);
auto &referenceE = inputWS->e(index_ref);
std::vector<double> refX(maxIt - minIt);
std::vector<double> refY(maxIt - minIt - 1);
std::vector<double> refE(maxIt - minIt - 1);
std::copy(minIt, maxIt, refX.begin());
mess << "min max" << refX.front() << " " << refX.back();
g_log.information(mess.str());
mess.str("");
std::copy(referenceY.cbegin() + difminIt,
referenceY.cbegin() + (difmaxIt - 1), refY.begin());
std::copy(referenceE.cbegin() + difminIt,
referenceE.cbegin() + (difmaxIt - 1), refE.begin());
// Now start the real stuff
// Create a 2DWorkspace that will hold the result
const int nY = static_cast<int>(refY.size());
const int npoints = 2 * nY - 3;
if (npoints < 1)
throw std::runtime_error("Range is not valid");
MatrixWorkspace_sptr out =
WorkspaceFactory::Instance().create(inputWS, nspecs, npoints, npoints);
// Calculate the mean value of the reference spectrum and associated error
// squared
double refMean = std::accumulate(refY.cbegin(), refY.cend(), 0.0);
double refMeanE2 = std::accumulate(refE.cbegin(), refE.cend(), 0.0,
VectorHelper::SumSquares<double>());
refMean /= static_cast<double>(nY);
refMeanE2 /= static_cast<double>(nY * nY);
auto itY = refY.begin();
auto itE = refE.begin();
double refVar = 0.0, refVarE = 0.0;
for (; itY != refY.end(); ++itY, ++itE) {
(*itY) -= refMean; // Now the vector is (y[i]-refMean)
(*itE) = (*itE) * (*itE) + refMeanE2; // New error squared
double t = (*itY) * (*itY); //(y[i]-refMean)^2
refVar += t; // Sum previous term
refVarE += 4.0 * t * (*itE); // Error squared
}
double refNorm = 1.0 / sqrt(refVar);
double refNormE = 0.5 * pow(refNorm, 3) * sqrt(refVarE);
// Now copy the other spectra
bool is_distrib = inputWS->isDistribution();
std::vector<double> XX(npoints);
for (int i = 0; i < npoints; ++i) {
XX[i] = static_cast<double>(i - nY + 2);
}
// Initialise the progress reporting object
out->mutableX(0) = XX;
m_progress = make_unique<Progress>(this, 0.0, 1.0, nspecs);
PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *out))
for (int i = 0; i < nspecs; ++i) // Now loop on all spectra
{
PARALLEL_START_INTERUPT_REGION
size_t wsIndex = indexes[i]; // Get the ws index from the table
// Copy spectra info from input Workspace
out->getSpectrum(i).copyInfoFrom(inputWS->getSpectrum(wsIndex));
out->setSharedX(i, out->sharedX(0));
// Get temp references
const auto &iX = inputWS->x(wsIndex);
const auto &iY = inputWS->y(wsIndex);
const auto &iE = inputWS->e(wsIndex);
// Copy Y,E data of spec(i) to temp vector
// Now rebin on the grid of reference spectrum
std::vector<double> tempY(nY);
std::vector<double> tempE(nY);
VectorHelper::rebin(iX.rawData(), iY.rawData(), iE.rawData(), refX, tempY,
tempE, is_distrib);
// Calculate the mean value of tempY
double tempMean = std::accumulate(tempY.begin(), tempY.end(), 0.0);
tempMean /= static_cast<double>(nY);
double tempMeanE2 = std::accumulate(tempE.begin(), tempE.end(), 0.0,
VectorHelper::SumSquares<double>());
tempMeanE2 /= static_cast<double>(nY * nY);
//
std::vector<double>::iterator itY;
std::vector<double>::iterator itE;
itY = tempY.begin();
itE = tempE.begin();
double tempVar = 0.0, tempVarE = 0.0;
for (; itY != tempY.end(); ++itY, ++itE) {
(*itY) -= tempMean; // Now the vector is (y[i]-refMean)
(*itE) = (*itE) * (*itE) + tempMeanE2; // New error squared double t = (*itY) * (*itY);
tempVar += t;
tempVarE += 4.0 * t * (*itE);
}
// Calculate the normalisation constant
double tempNorm = 1.0 / sqrt(tempVar);
double tempNormE = 0.5 * pow(tempNorm, 3) * sqrt(tempVarE);
double normalisation = refNorm * tempNorm;
double normalisationE2 =
pow((refNorm * tempNormE), 2) + pow((tempNorm * refNormE), 2);
// Get reference to the ouput spectrum
auto &outY = out->mutableY(i);
auto &outE = out->mutableE(i);
for (int k = -nY + 2; k <= nY - 2; ++k) {
int kp = abs(k);
double val = 0, err2 = 0, x, y, xE, yE;
for (int j = nY - 1 - kp; j >= 0; --j) {
if (k >= 0) {
x = refY[j];
y = tempY[j + kp];
xE = refE[j];
yE = tempE[j + kp];
} else {
x = tempY[j];
y = refY[j + kp];
xE = tempE[j];
yE = refE[j + kp];
}
val += (x * y);
err2 += x * x * yE + y * y * xE;
}
outY[k + nY - 2] = (val * normalisation);
outE[k + nY - 2] = sqrt(val * val * normalisationE2 +
normalisation * normalisation * err2);
}
// Update progress information
// double prog=static_cast<double>(i)/nspecs;
// progress(prog);
m_progress->report();
// interruption_point();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
setProperty("OutputWorkspace", out);
}
} // namespace Algorithm
} // namespace Mantid