Commit d90cb366 authored by Zhou, Wenduo's avatar Zhou, Wenduo
Browse files

Implement rwp calculator. Refs #6746.

parent 6bdca971
......@@ -62,6 +62,8 @@ namespace Mantid
/// Return the modified Z score values for a dataset
template<typename TYPE>
std::vector<double> getModifiedZscore(const std::vector<TYPE>& data, const bool sorted=false);
/// Return the R-factors (Rwp) of a diffraction pattern data
double getRFactor(const std::vector<double>& obsI, const std::vector<double>& calI, const std::vector<double>& obsE);
} // namespace Kernel
} // namespace Mantid
......
......@@ -7,6 +7,8 @@
#include <cmath>
#include <numeric>
#include <string>
#include <stdexcept>
#include <sstream>
namespace Mantid
{
......@@ -214,6 +216,43 @@ namespace Mantid
return getNanStatistics();
}
/** Return the Rwp of a diffraction pattern data
* @param obsY :: array of observed intensity values
* @param calY :: array of calculated intensity values;
*
*/
double getRFactor(const std::vector<double>& obsI, const std::vector<double>& calI, const std::vector<double>& obsE)
{
// 1. Check
if (obsI.size() != calI.size() || obsI.size() != obsE.size())
{
std::stringstream errss;
errss << "GetRFactor() Input Error! Observed Intensity, Calculated Intensity and Observed Error have different number of elements.";
throw std::runtime_error(errss.str());
}
double sumnom = 0;
double sumdenom = 0;
size_t numpts = obsI.size();
for (size_t i = 0; i < numpts; ++i)
{
double cal_i = calI[i];
double obs_i = obsI[i];
double sigma = obsE[i];
double weight = 1.0/(sigma*sigma);
double diff = obs_i - cal_i;
sumnom += weight*diff*diff;
sumdenom += weight*obs_i*obs_i;
}
double rwp = std::sqrt(sumnom/sumdenom);
return rwp;
}
// -------------------------- Macro to instantiation concrete types --------------------------------
#define INSTANTIATE(TYPE) \
template MANTID_KERNEL_DLL Statistics getStatistics<TYPE> (const vector<TYPE> &, const bool); \
......
......@@ -115,6 +115,75 @@ public:
TS_ASSERT(my_isnan(stats.maximum));
TS_ASSERT(my_isnan(stats.median));
}
/** Test function to calculate Rwp
*/
void testRwp()
{
vector<double> obsY(4);
vector<double> calY(4);
vector<double> obsE(4);
obsY[0] = 1.0;
calY[0] = 1.1;
obsE[0] = 1.0;
obsY[1] = 2.0;
calY[1] = 2.1;
obsE[1] = 1.2;
obsY[2] = 3.0;
calY[2] = 3.5;
obsE[2] = 1.4;
obsY[3] = 1.0;
calY[3] = 1.3;
obsE[3] = 1.0;
double rwp = getRFactor(obsY, calY, obsE);
TS_ASSERT_DELTA(rwp, 0.1582, 0.0001);
}
};
#endif // STATISTICSTEST_H_
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment