14#include "MantidHistogramData/Histogram.h"
15#include "MantidHistogramData/LinearGenerator.h"
31using namespace HistogramData;
38using namespace DataObjects;
42const string BIG_G_OF_R(
"G(r)");
44const string LITTLE_G_OF_R(
"g(r)");
46const string RDF_OF_R(
"RDF(r)");
48const string G_K_OF_R(
"G_k(r)");
51const string S_OF_Q(
"S(Q)");
53const string S_OF_Q_MINUS_ONE(
"S(Q)-1");
55const string Q_S_OF_Q_MINUS_ONE(
"Q[S(Q)-1]");
57const string FORWARD(
"Forward");
59const string BACKWARD(
"Backward");
72 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
74 "Input spectrum density or paired-distribution function");
76 "Result paired-distribution function or Input spectrum density");
78 std::vector<std::string> directionOptions;
79 directionOptions.emplace_back(FORWARD);
80 directionOptions.emplace_back(BACKWARD);
81 declareProperty(
"Direction", FORWARD, std::make_shared<StringListValidator>(directionOptions),
82 "The direction of the fourier transform");
84 "Average number density used for g(r) and RDF(r) conversions (optional)");
85 declareProperty(
"Filter",
false,
"Set to apply Lorch function filter to the input");
88 std::vector<std::string> soqTypes;
89 soqTypes.emplace_back(S_OF_Q);
90 soqTypes.emplace_back(S_OF_Q_MINUS_ONE);
91 soqTypes.emplace_back(Q_S_OF_Q_MINUS_ONE);
92 declareProperty(
"InputSofQType", S_OF_Q, std::make_shared<StringListValidator>(soqTypes),
93 "To identify spectral density function (deprecated)");
95 declareProperty(
"SofQType", S_OF_Q, std::make_shared<StringListValidator>(soqTypes),
96 "To identify spectral density function");
97 mustBePositive->setLower(0.0);
100 "Step size of Q of S(Q) to calculate. Default = "
101 ":math:`\\frac{\\pi}{R_{max}}`.");
104 "Minimum Q in S(Q) to calculate in Fourier transform (optional).");
107 "Maximum Q in S(Q) to calculate in Fourier transform. "
108 "(optional, defaults to 40 on backward transform.)");
111 std::vector<std::string> pdfTypes;
112 pdfTypes.emplace_back(BIG_G_OF_R);
113 pdfTypes.emplace_back(LITTLE_G_OF_R);
114 pdfTypes.emplace_back(RDF_OF_R);
115 pdfTypes.emplace_back(G_K_OF_R);
116 declareProperty(
"PDFType", BIG_G_OF_R, std::make_shared<StringListValidator>(pdfTypes),
117 "Type of output PDF including G(r)");
120 "Step size of r of G(r) to calculate. Default = "
121 ":math:`\\frac{\\pi}{Q_{max}}`.");
126 "Maximum r for G(r) to calculate. (optional, defaults to 20 "
127 "on forward transform.)");
129 string recipGroup(
"Reciprocal Space");
135 string realGroup(
"Real Space");
143 std::map<string, string> result;
149 result[
"Qmax"] =
"Must be greater than Qmin";
157 if (inputWS->getNumberHistograms() != 1) {
158 result[
"InputWorkspace"] =
"Input workspace must have only one spectrum";
160 const std::string inputXunit = inputWS->getAxis(0)->unit()->unitID();
161 if (inputXunit !=
"MomentumTransfer" && inputXunit !=
"dSpacing" && inputXunit !=
"AtomicDistance") {
162 result[
"InputWorkspace"] =
"Input workspace units not supported";
172 }
else if (min <
X.front()) {
173 g_log.
information(
"Specified input min < range of data. Adjusting to data range.");
178 auto iter = std::lower_bound(
X.begin(),
X.end(), min);
179 size_t min_index = std::distance(
X.begin(), iter);
182 auto iter_first_normal =
183 std::find_if(std::next(
Y.begin(), min_index),
Y.end(),
static_cast<bool (*)(
double)
>(std::isnormal));
184 size_t first_normal_index = std::distance(
Y.begin(), iter_first_normal);
185 if (first_normal_index > min_index) {
186 g_log.
information(
"Specified input min where data is nan/inf or zero. Adjusting to data range.");
187 min_index = first_normal_index;
197 }
else if (max >
X.back()) {
198 g_log.
information() <<
"Specified input max > range of data. Adjusting to data range.\n";
203 auto iter = std::upper_bound(
X.begin(),
X.end(), max) - 1;
204 size_t max_index = std::distance(
X.begin(), iter);
209 auto back_iter = std::find_if(
Y.rbegin(),
Y.rend(),
static_cast<bool (*)(
double)
>(std::isnormal));
210 size_t first_normal_index =
Y.size() - std::distance(
Y.rbegin(), back_iter);
211 if (first_normal_index < max_index) {
212 g_log.
information(
"Specified input max where data is nan/inf or zero. Adjusting to data range.");
213 max_index = first_normal_index;
228 if (!
isEmpty(materialDensity) && materialDensity > 0)
229 rho0 = materialDensity;
238 std::vector<double> &DFOfQ,
const std::vector<double> &DQ) {
241 string inputSOQType =
getProperty(
"InputSofQType");
243 soqType = inputSOQType;
244 g_log.
warning() <<
"InputSofQType has been deprecated and replaced by SofQType\n";
246 if (soqType == S_OF_Q) {
249 std::for_each(FOfQ.begin(), FOfQ.end(), [](
double &F) { F--; });
250 soqType = S_OF_Q_MINUS_ONE;
252 if (soqType == Q_S_OF_Q_MINUS_ONE) {
255 for (
size_t i = 0; i < DFOfQ.size(); ++i) {
256 DFOfQ[i] = (Q[i] / DQ[i] + FOfQ[i] / DFOfQ[i]) * (FOfQ[i] / Q[i]);
259 std::transform(FOfQ.begin(), FOfQ.end(), FOfQ.begin(), Q.begin(), std::divides<double>());
260 soqType = S_OF_Q_MINUS_ONE;
262 if (soqType != S_OF_Q_MINUS_ONE) {
264 std::stringstream msg;
265 msg <<
"Do not understand SofQType = " << soqType;
266 throw std::runtime_error(msg.str());
272 std::vector<double> &DFOfR,
const std::vector<double> &DR,
273 const std::string &PDFType,
const double &rho0,
274 const double &cohScatLen) {
275 if (PDFType == LITTLE_G_OF_R) {
276 for (
size_t i = 0; i < FOfR.size(); ++i) {
278 FOfR[i] = FOfR[i] - 1.0;
280 }
else if (PDFType == BIG_G_OF_R) {
281 const double factor = 4. * M_PI * rho0;
282 for (
size_t i = 0; i < FOfR.size(); ++i) {
284 DFOfR[i] = (R[i] / DR[i] + FOfR[i] / DFOfR[i]) * (FOfR[i] / R[i]);
286 FOfR[i] = FOfR[i] / (factor * R[i]);
288 }
else if (PDFType == RDF_OF_R) {
289 const double factor = 4. * M_PI * rho0;
290 for (
size_t i = 0; i < FOfR.size(); ++i) {
292 DFOfR[i] = (2.0 * R[i] / DR[i] + FOfR[i] / DFOfR[i]) * (FOfR[i] / R[i]);
294 FOfR[i] = FOfR[i] / (factor * R[i] * R[i]) - 1.0;
296 }
else if (PDFType == G_K_OF_R) {
297 const double factor = 0.01 * pow(cohScatLen, 2);
299 for (
size_t i = 0; i < FOfR.size(); ++i) {
301 DFOfR[i] = DFOfR[i] / factor;
303 FOfR[i] = FOfR[i] / factor;
310 HistogramData::HistogramE &DFOfQ) {
313 string inputSOQType =
getProperty(
"InputSofQType");
315 soqType = inputSOQType;
316 g_log.
warning() <<
"InputSofQType has been deprecated and replaced by SofQType\n";
318 if (soqType == S_OF_Q) {
319 for (
size_t i = 0; i < FOfQ.size(); ++i) {
321 FOfQ[i] = FOfQ[i] + 1.0;
323 }
else if (soqType == Q_S_OF_Q_MINUS_ONE) {
324 for (
size_t i = 0; i < FOfQ.size(); ++i) {
325 DFOfQ[i] = Q[i] * DFOfQ[i];
326 FOfQ[i] = FOfQ[i] * Q[i];
333 const HistogramData::HistogramX &R,
334 HistogramData::HistogramE &DFOfR,
const std::string &PDFType,
335 const double &rho0,
const double &cohScatLen) {
337 if (PDFType == LITTLE_G_OF_R) {
338 for (
size_t i = 0; i < FOfR.size(); ++i) {
340 FOfR[i] = FOfR[i] + 1.0;
342 }
else if (PDFType == BIG_G_OF_R) {
343 const double factor = 4. * M_PI * rho0;
344 for (
size_t i = 0; i < FOfR.size(); ++i) {
346 DFOfR[i] = DFOfR[i] * R[i];
348 FOfR[i] = FOfR[i] * factor * R[i];
350 }
else if (PDFType == RDF_OF_R) {
351 const double factor = 4. * M_PI * rho0;
352 for (
size_t i = 0; i < FOfR.size(); ++i) {
354 DFOfR[i] = DFOfR[i] * R[i];
356 FOfR[i] = (FOfR[i] + 1.0) * factor * R[i] * R[i];
358 }
else if (PDFType == G_K_OF_R) {
359 const double factor = 0.01 * pow(cohScatLen, 2);
361 for (
size_t i = 0; i < FOfR.size(); ++i) {
363 DFOfR[i] = DFOfR[i] * factor;
365 FOfR[i] = FOfR[i] * factor;
376 auto inputX = inputWS->x(0).rawData();
377 std::vector<double> inputDX(inputX.size(), 0.0);
378 if (inputWS->sharedDx(0))
379 inputDX = inputWS->dx(0).rawData();
380 auto inputY = inputWS->y(0).rawData();
381 auto inputDY = inputWS->e(0).rawData();
385 const std::string inputXunit = inputWS->getAxis(0)->unit()->unitID();
386 if (inputXunit ==
"MomentumTransfer") {
388 }
else if (inputXunit ==
"dSpacing") {
390 const double PI_2(2. * M_PI);
391 std::for_each(inputX.begin(), inputX.end(), [&PI_2](
double &Q) { Q /= PI_2; });
392 std::transform(inputDX.begin(), inputDX.end(), inputX.begin(), inputDX.begin(), std::divides<double>());
394 std::reverse(inputX.begin(), inputX.end());
395 std::reverse(inputDX.begin(), inputDX.end());
396 std::reverse(inputY.begin(), inputY.end());
397 std::reverse(inputDY.begin(), inputDY.end());
398 }
else if (inputXunit ==
"AtomicDistance") {
401 g_log.
debug() <<
"Input unit is " << inputXunit <<
"\n";
404 if (!inputWS->isHistogramData()) {
405 g_log.
warning() <<
"This algorithm has not been tested on density data "
406 "(only on histograms)\n";
422 const std::string PDFType =
getProperty(
"PDFType");
428 if (PDFType ==
"G_k(r)" && cohScatLen == 0.0) {
429 std::stringstream msg;
430 msg <<
"Coherent Scattering Length is zero. Please check a sample material has been specified";
431 throw std::runtime_error(msg.str());
435 if (direction == FORWARD) {
437 }
else if (direction == BACKWARD) {
441 double inMin, inMax, outDelta, outMax;
449 if (direction == BACKWARD) {
462 g_log.
notice() <<
"Adjusting to data: input min = " << inputX[Xmin_index] <<
" input max = " << inputX[Xmax_index]
466 outDelta = M_PI / inputX[Xmax_index];
467 auto sizer =
static_cast<size_t>(outMax / outDelta);
474 outputWS->copyExperimentInfoFrom(inputWS.get());
475 if (direction == FORWARD) {
477 outputWS->setYUnitLabel(
"PDF");
478 outputWS->mutableRun().addProperty(
"Qmin", inputX[Xmin_index],
"Angstroms^-1",
true);
479 outputWS->mutableRun().addProperty(
"Qmax", inputX[Xmax_index],
"Angstroms^-1",
true);
480 }
else if (direction == BACKWARD) {
482 outputWS->setYUnitLabel(
"Spectrum Density");
483 outputWS->mutableRun().addProperty(
"Rmin", inputX[Xmin_index],
"Angstroms",
true);
484 outputWS->mutableRun().addProperty(
"Rmax", inputX[Xmax_index],
"Angstroms",
true);
486 outputWS->setDistribution(
true);
487 BinEdges edges(sizer + 1, LinearGenerator(outDelta, outDelta));
488 outputWS->setBinEdges(0, edges);
489 auto &outputX = outputWS->mutableX(0);
490 g_log.
information() <<
"Using output min = " << outputX.front() <<
"and output max = " << outputX.back() <<
"\n";
492 auto &outputY = outputWS->mutableY(0);
493 auto &outputE = outputWS->mutableE(0);
499 double corr = 0.5 / M_PI / M_PI / rho0;
500 if (direction == BACKWARD) {
501 corr = 4.0 * M_PI * rho0;
503 for (
size_t outXIndex = 0; outXIndex < sizer; outXIndex++) {
504 const double outX = outputX[outXIndex];
505 const double outXFac = corr / (outX * outX * outX);
508 double errorSquared = 0;
509 auto pdfIntegral = [outX](
const double x) ->
double {
return sin(
x * outX) -
x * outX * cos(
x * outX); };
510 double inX2 = inputX[Xmin_index];
511 double integralX2 = pdfIntegral(inX2);
512 const double inXMax = inputX[Xmax_index];
514 for (
size_t inXIndex = Xmin_index; inXIndex < Xmax_index; inXIndex++) {
515 const double inX1 = inX2;
516 inX2 = inputX[inXIndex + 1];
517 const double integralX1 = integralX2;
518 integralX2 = pdfIntegral(inX2);
519 double defIntegral = integralX2 - integralX1;
522 if (filter && inX1 != 0) {
523 const double lorchKernel = inX1 * M_PI / inXMax;
524 defIntegral *= sin(lorchKernel) / lorchKernel;
526 fs += defIntegral * inputY[inXIndex];
528 const double error = defIntegral * inputDY[inXIndex];
533 outputY[outXIndex] = fs * outXFac;
534 outputE[outXIndex] = sqrt(errorSquared) * outXFac;
537 if (direction == FORWARD) {
539 }
else if (direction == BACKWARD) {
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
bool isDefault(const std::string &name) const
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void debug(const std::string &msg)
Logs at debug level.
void notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAt...
double cohScatterLength(const double lambda=PhysicalConstants::NeutronAtom::ReferenceLambda) const
Get the coherent scattering length for a given wavelength in fm.
double numberDensity() const
Get the number density.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.