Skip to content
Snippets Groups Projects
Commit df894b29 authored by Anton Piccardo-Selg's avatar Anton Piccardo-Selg
Browse files

Refs #11534 Implement dx for Q1D

parent b1f26b42
No related branches found
No related tags found
No related merge requests found
...@@ -42,6 +42,10 @@ public: ...@@ -42,6 +42,10 @@ public:
API::MatrixWorkspace_const_sptr detectAdj, API::MatrixWorkspace_const_sptr detectAdj,
API::MatrixWorkspace_const_sptr qResolution); API::MatrixWorkspace_const_sptr qResolution);
void examineInput(API::MatrixWorkspace_const_sptr dataWS,
API::MatrixWorkspace_const_sptr binAdj,
API::MatrixWorkspace_const_sptr detectAdj);
size_t waveLengthCutOff(API::MatrixWorkspace_const_sptr dataWS, size_t waveLengthCutOff(API::MatrixWorkspace_const_sptr dataWS,
const double RCut, const double WCut, const double RCut, const double WCut,
const size_t specInd) const; const size_t specInd) const;
......
...@@ -129,6 +129,8 @@ void Q1D2::exec() { ...@@ -129,6 +129,8 @@ void Q1D2::exec() {
MantidVec normSum(YOut.size(), 0.0); MantidVec normSum(YOut.size(), 0.0);
// the error on the normalisation // the error on the normalisation
MantidVec normError2(YOut.size(), 0.0); MantidVec normError2(YOut.size(), 0.0);
// the averaged Q resolution
MantidVec &qResolutionOut = outputWS->dataDx(0);
const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms()); const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
Progress progress(this, 0.05, 1.0, numSpec + 1); Progress progress(this, 0.05, 1.0, numSpec + 1);
...@@ -167,6 +169,7 @@ void Q1D2::exec() { ...@@ -167,6 +169,7 @@ void Q1D2::exec() {
continue; continue;
} }
const size_t numWavbins = m_dataWS->readY(i).size() - wavStart; const size_t numWavbins = m_dataWS->readY(i).size() - wavStart;
// make just one call to new to reduce CPU overhead on each thread, access // make just one call to new to reduce CPU overhead on each thread, access
// to these // to these
...@@ -191,13 +194,17 @@ void Q1D2::exec() { ...@@ -191,13 +194,17 @@ void Q1D2::exec() {
MantidVec::const_iterator YIn = m_dataWS->readY(i).begin() + wavStart; MantidVec::const_iterator YIn = m_dataWS->readY(i).begin() + wavStart;
MantidVec::const_iterator EIn = m_dataWS->readE(i).begin() + wavStart; MantidVec::const_iterator EIn = m_dataWS->readE(i).begin() + wavStart;
// Pointers to the QResolution data. Note that the xdata was initially the same, hence
// the same indexing applies to the y values of m_dataWS and qResolution
MantidVec::const_iterator QResIn = qResolution->readY(i).begin() + wavStart;
// when finding the output Q bin remember that the input Q bins (from the // when finding the output Q bin remember that the input Q bins (from the
// convert to wavelength) start high and reduce // convert to wavelength) start high and reduce
MantidVec::const_iterator loc = QOut.end(); MantidVec::const_iterator loc = QOut.end();
// sum the Q contributions from each individual spectrum into the output // sum the Q contributions from each individual spectrum into the output
// array // array
const MantidVec::const_iterator end = m_dataWS->readY(i).end(); const MantidVec::const_iterator end = m_dataWS->readY(i).end();
for (; YIn != end; ++YIn, ++EIn, ++QIn, ++norms, ++normETo2s) { for (; YIn != end; ++YIn, ++QResIn, ++EIn, ++QIn, ++norms, ++normETo2s) {
// find the output bin that each input y-value will fall into, remembering // find the output bin that each input y-value will fall into, remembering
// there is one more bin boundary than bins // there is one more bin boundary than bins
getQBinPlus1(QOut, *QIn, loc); getQBinPlus1(QOut, *QIn, loc);
...@@ -212,13 +219,13 @@ void Q1D2::exec() { ...@@ -212,13 +219,13 @@ void Q1D2::exec() {
// at the end // at the end
EOutTo2[bin] += (*EIn) * (*EIn); EOutTo2[bin] += (*EIn) * (*EIn);
normError2[bin] += *normETo2s; normError2[bin] += *normETo2s;
qResolutionOut[bin] += (*YIn) *(*QResIn);
} }
} }
} }
PARALLEL_CRITICAL(q1d_spectra_map) { PARALLEL_CRITICAL(q1d_spectra_map) {
progress.report("Computing I(Q)"); progress.report("Computing I(Q)");
// Add up the detector IDs in the output spectrum at workspace index 0 // Add up the detector IDs in the output spectrum at workspace index 0
const ISpectrum *inSpec = m_dataWS->getSpectrum(i); const ISpectrum *inSpec = m_dataWS->getSpectrum(i);
ISpectrum *outSpec = outputWS->getSpectrum(0); ISpectrum *outSpec = outputWS->getSpectrum(0);
...@@ -229,12 +236,21 @@ void Q1D2::exec() { ...@@ -229,12 +236,21 @@ void Q1D2::exec() {
} }
PARALLEL_CHECK_INTERUPT_REGION PARALLEL_CHECK_INTERUPT_REGION
// Calculate the QResolution and add it as a DX value to the outputworkspace
Mantid::MantidVec::const_iterator countsIterator = YOut.begin();
Mantid::MantidVec::iterator qResolutionIterator = qResolutionOut.begin();
for (;countsIterator != YOut.end(); ++countsIterator, ++qResolutionIterator) {
// Divide by the counts of the Qbin
*qResolutionIterator = (*qResolutionIterator)/(*countsIterator);
}
bool doOutputParts = getProperty("OutputParts"); bool doOutputParts = getProperty("OutputParts");
if (doOutputParts) { if (doOutputParts) {
MatrixWorkspace_sptr ws_sumOfCounts = MatrixWorkspace_sptr ws_sumOfCounts =
WorkspaceFactory::Instance().create(outputWS); WorkspaceFactory::Instance().create(outputWS);
ws_sumOfCounts->dataX(0) = outputWS->dataX(0); ws_sumOfCounts->dataX(0) = outputWS->dataX(0);
ws_sumOfCounts->dataY(0) = outputWS->dataY(0); ws_sumOfCounts->dataY(0) = outputWS->dataY(0);
ws_sumOfCounts->dataDx(0) = outputWS->dataDx(0);
for (size_t i = 0; i < outputWS->dataE(0).size(); i++) { for (size_t i = 0; i < outputWS->dataE(0).size(); i++) {
ws_sumOfCounts->dataE(0)[i] = sqrt(outputWS->dataE(0)[i]); ws_sumOfCounts->dataE(0)[i] = sqrt(outputWS->dataE(0)[i]);
} }
...@@ -242,6 +258,7 @@ void Q1D2::exec() { ...@@ -242,6 +258,7 @@ void Q1D2::exec() {
MatrixWorkspace_sptr ws_sumOfNormFactors = MatrixWorkspace_sptr ws_sumOfNormFactors =
WorkspaceFactory::Instance().create(outputWS); WorkspaceFactory::Instance().create(outputWS);
ws_sumOfNormFactors->dataX(0) = outputWS->dataX(0); ws_sumOfNormFactors->dataX(0) = outputWS->dataX(0);
ws_sumOfNormFactors->dataDx(0) = outputWS->dataDx(0);
for (size_t i = 0; i < ws_sumOfNormFactors->dataY(0).size(); i++) { for (size_t i = 0; i < ws_sumOfNormFactors->dataY(0).size(); i++) {
ws_sumOfNormFactors->dataY(0)[i] = normSum[i]; ws_sumOfNormFactors->dataY(0)[i] = normSum[i];
ws_sumOfNormFactors->dataE(0)[i] = sqrt(normError2[i]); ws_sumOfNormFactors->dataE(0)[i] = sqrt(normError2[i]);
......
...@@ -31,6 +31,44 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS, ...@@ -31,6 +31,44 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
API::MatrixWorkspace_const_sptr binAdj, API::MatrixWorkspace_const_sptr binAdj,
API::MatrixWorkspace_const_sptr detectAdj, API::MatrixWorkspace_const_sptr detectAdj,
API::MatrixWorkspace_const_sptr qResolution) { API::MatrixWorkspace_const_sptr qResolution) {
// Check the compatibility of dataWS, binAdj and detectAdj
examineInput(dataWS, binAdj, detectAdj);
// Check the compatibility of the QResolution workspace
if (qResolution) {
// We require the same number of histograms
if (detectAdj->getNumberHistograms() != dataWS->getNumberHistograms()) {
throw std::invalid_argument("The QResolution should have one spectrum"
"per spectrum of the input workspace");
}
// We require the same binning for the input workspace and the q resolution
// workspace
MantidVec::const_iterator reqX = dataWS->readX(0).begin();
MantidVec::const_iterator qResX = qResolution->readX(0).begin();
for (; reqX != dataWS->readX(0).end(); ++reqX, ++qResX) {
if (*reqX != *qResX) {
throw std::invalid_argument("The QResolution needs to have the same binning as"
"as the input workspace.");
}
}
}
}
/** Checks if workspaces input to Q1D or Qxy are reasonable
@param dataWS data workspace
@param binAdj (WavelengthAdj) workpace that will be checked to see if it has
one spectrum and the same number of bins as dataWS
@param detectAdj (PixelAdj) passing NULL for this wont raise an error, if set
it will be checked this workspace has as many histograms as dataWS each with
one bin
@throw invalid_argument if the workspaces are not mututially compatible
*/
void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
API::MatrixWorkspace_const_sptr binAdj,
API::MatrixWorkspace_const_sptr detectAdj) {
if (dataWS->getNumberHistograms() < 1) { if (dataWS->getNumberHistograms() < 1) {
throw std::invalid_argument( throw std::invalid_argument(
"Empty data workspace passed, can not continue"); "Empty data workspace passed, can not continue");
...@@ -103,28 +141,9 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS, ...@@ -103,28 +141,9 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
} }
} }
} }
// Check the compatibility of the QResolution workspace
if (qResolution) {
// We require the same number of histograms
if (detectAdj->getNumberHistograms() != dataWS->getNumberHistograms()) {
throw std::invalid_argument("The QResolution should have one spectrum"
"per spectrum of the input workspace");
}
// We require the same binning for the input workspace and the q resolution
// workspace
MantidVec::const_iterator reqX = dataWS->readX(0).begin();
MantidVec::const_iterator qResX = qResolution->readX(0).begin();
for (; reqX != dataWS->readX(0).end(); ++reqX, ++qResX) {
if (*reqX != *qResX) {
throw std::invalid_argument("The QResolution needs to have the same binning as"
"as the input workspace.");
}
}
}
} }
/** Finds the first index number of the first wavelength bin that should /** Finds the first index number of the first wavelength bin that should
* included based on the * included based on the
* the calculation: W = Wcut (Rcut-R)/Rcut * the calculation: W = Wcut (Rcut-R)/Rcut
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment