Skip to content

Commit

Permalink
Started to refactor the unit conversion path. Refs #10929.
Browse files Browse the repository at this point in the history
  • Loading branch information
wdzhou committed Feb 21, 2015
1 parent 68526f9 commit a6ab212
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 24 deletions.
Expand Up @@ -16,8 +16,22 @@ namespace MDAlgorithms {
* @param wavelength
* @return
*/
inline double calculateD(const double &theta, const double &wavelength) {
return (0.5 * wavelength / sin(theta));
inline double calculateDspaceFrom2Theta(const double &twotheta,
const double &wavelength) {
return (0.5 * wavelength / sin(twotheta * 0.5 * M_PI / 180.));
}

//----------------------------------------------------------------------------------------------
/** Calcualte detector's 2theta value from d-spacing value
* 2theta = 2*asin(lamba/2/d)
* @brief calculate2ThetaFromD
* @param dspace
* @param wavelengh
* @return
*/
inline double calculate2ThetaFromD(const double &dspace,
const double &wavelength) {
return (2.0 * asin(0.5 * wavelength / dspace) * 180. / M_PI);
}

//----------------------------------------------------------------------------------------------
Expand All @@ -28,8 +42,21 @@ inline double calculateD(const double &theta, const double &wavelength) {
* @param wavelength
* @return
*/
inline double calculateQ(const double &theta, const double &wavelength) {
return (4 * M_PI * sin(theta) / wavelength);
inline double calculateQFrom2Theta(const double &twotheta,
const double &wavelength) {
return (4. * M_PI * sin(twotheta * 0.5 * M_PI / 180.) / wavelength);
}

//----------------------------------------------------------------------------------------------
/** Calculate detector's 2theta value from Q
* 2theta = 2*asin(lambda*Q/(4pi))
* @brief calcualte2ThetaFromQ
* @param q
* @param wavelength
* @return
*/
inline double calcualte2ThetaFromQ(const double &q, const double &wavelength) {
return (2.0 * asin(q * wavelength * 0.25 / M_PI) * 180. / M_PI);
}

/** ConvertCWPDMDToSpectra : Convert one MDWorkspaces containing reactor-source
Expand Down
83 changes: 63 additions & 20 deletions Code/Mantid/Framework/MDAlgorithms/src/ConvertCWPDMDToSpectra.cpp
Expand Up @@ -95,14 +95,74 @@ void ConvertCWPDMDToSpectra::exec() {
double scaleFactor = getProperty("ScaleFactor");
// do linear interpolation
bool doLinearInterpolation = getProperty("LinearInterpolateZeroCounts");
// unit
std::string outputunit = getProperty("UnitOutput");
double wavelength = getProperty("NeutronWaveLength");

// Validate inputs
// input data workspace and monitor workspace should match
size_t numdataevents = inputDataWS->getNEvents();
size_t nummonitorevents = inputMonitorWS->getNEvents();
if (numdataevents != nummonitorevents)
throw std::runtime_error("Input data workspace and monitor workspace have "
"different number of MDEvents.");

// output unit: make a map for wavelength
std::map<uint16_t, double> map_runWavelength;
if (outputunit.compare("2theta")) {
// set up runid and wavelength map
std::string wavelengthpropertyname =
getProperty("NeutornWaveLengthPropertyName");

uint16_t numexpinfo = inputDataWS->getNumExperimentInfo();
for (uint16_t iexp = 0; iexp < numexpinfo; ++iexp) {
int runid = atoi(inputDataWS->getExperimentInfo(iexp)
->run()
.getProperty("run")
->value()
.c_str());
double thislambda;
if (inputDataWS->getExperimentInfo(iexp)->run().hasProperty(
wavelengthpropertyname))
thislambda = atof(inputDataWS->getExperimentInfo(iexp)
->run()
.getProperty(wavelengthpropertyname)
->value()
.c_str());
else if (wavelength != EMPTY_DBL())
thislambda = wavelength;
else {
std::stringstream errss;
errss << "In order to convert unit to " << outputunit
<< ", either NeutronWaveLength "
" is to be specified or property " << wavelengthpropertyname
<< " must exist for run " << runid << ".";
throw std::runtime_error(errss.str());
}
}
}

/*
if ( outputunit.compare("2theta") && (wavelength == EMPTY_DBL()) ) {
// Find it via property
std::string wavelengthpropertyname =
getProperty("NeutornWaveLengthPropertyName");
if
(inputDataWS->getExperimentInfo(0)->run().hasProperty(wavelengthpropertyname))
{
std::stringstream errss;
errss << "In order to convert unit to " << outputunit
<< ", either NeutronWaveLength "
" is to be specified or property " << wavelengthpropertyname
<< " must exist.";
throw std::runtime_error(errss.str());
}
wavelength = atof(
outws->run().getProperty(wavelengthpropertyname)->value().c_str());
}
*/

// bin parameters
if (binParams.size() != 3)
throw std::runtime_error("Binning parameters must have 3 items.");
if (binParams[0] >= binParams[2])
Expand All @@ -124,24 +184,7 @@ void ConvertCWPDMDToSpectra::exec() {
<< "\n";

// Output units
std::string outputunit = getProperty("UnitOutput");
if (outputunit.compare("2theta")) {
double wavelength = getProperty("NeutronWaveLength");
if (wavelength == EMPTY_DBL()) {
// Find it via property
std::string wavelengthpropertyname =
getProperty("NeutornWaveLengthPropertyName");
if (!outws->run().hasProperty(wavelengthpropertyname)) {
std::stringstream errss;
errss << "In order to convert unit to " << outputunit
<< ", either NeutronWaveLength "
" is to be specified or property " << wavelengthpropertyname
<< " must exist.";
throw std::runtime_error(errss.str());
}
wavelength = atof(
outws->run().getProperty(wavelengthpropertyname)->value().c_str());
}
convertUnits(outws, outputunit, wavelength);
}
// TODO : Implement unit for 2theta in another ticket.
Expand Down Expand Up @@ -312,7 +355,7 @@ void ConvertCWPDMDToSpectra::binMD(API::IMDEventWorkspace_const_sptr mdws,
if (xindex < 0)
g_log.warning("xindex < 0");
if (xindex >= static_cast<int>(vecy.size()) - 1) {
// FIXME - It may throw away all the detectors' value above Ymax
// If the Xmax is set too narrow, then
g_log.error() << "This is the bug! "
<< "xindex = " << xindex << " 2theta = " << twotheta
<< " out of [" << vecx.front() << ", " << vecx.back()
Expand Down Expand Up @@ -494,9 +537,9 @@ void ConvertCWPDMDToSpectra::convertUnits(API::MatrixWorkspace_sptr matrixws,
size_t vecsize = vecX.size();
for (size_t j = 0; j < vecsize; ++j) {
if (target == 'd')
vecX[j] = calculateD(vecX[j] * 0.5, wavelength);
vecX[j] = calculateDspaceFrom2Theta(vecX[j] * 0.5, wavelength);
else
vecX[j] = calculateQ(vecX[j] * 0.5, wavelength);
vecX[j] = calculateQFrom2Theta(vecX[j] * 0.5, wavelength);
}
}

Expand Down

0 comments on commit a6ab212

Please sign in to comment.