From 1eada25135e90cde37b864edf6335e6a90dcaf6b Mon Sep 17 00:00:00 2001
From: Steve Williams <stephen.williams@stfc.ac.uk>
Date: Fri, 24 Jun 2011 19:35:41 +0000
Subject: [PATCH] Q1D has been updated! re #3167 fixes #3105
---
.../Algorithms/inc/MantidAlgorithms/Q1D2.h | 19 +-
Code/Mantid/Framework/Algorithms/src/Q1D2.cpp | 341 ++++++++++++------
.../Framework/Algorithms/test/Q1D2Test.h | 44 ++-
3 files changed, 290 insertions(+), 114 deletions(-)
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Q1D2.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Q1D2.h
index 5b012dc6e86..cdaace0d6c3 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Q1D2.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Q1D2.h
@@ -41,7 +41,7 @@ class DLLExport Q1D2 : public API::Algorithm
{
public:
/// (Empty) Constructor
- Q1D2() : API::Algorithm(), m_numInBins(0) {}
+ Q1D2() : API::Algorithm(), m_RCut(0.0), m_WCutOver(0.0) {}
/// Virtual destructor
virtual ~Q1D2() {}
/// Algorithm's name
@@ -54,8 +54,10 @@ public:
private:
/// the experimental workspace with counts across the detector
API::MatrixWorkspace_const_sptr m_dataWS;
- /// number of bin boundies in the input workspace
- size_t m_numInBins;
+ ///The radius cut off, should be value of the property RadiusCut. A value of zero here will disable the cut off and all wavelengths are used
+ double m_RCut;
+ ///The wavelength cut off divided by the radius cut is used in the calculation of the first wavelength to include, it's value is only used if RadiusCut > 0
+ double m_WCutOver;
/// Sets documentation strings for this algorithm
virtual void initDocs();
@@ -64,14 +66,19 @@ private:
/// Execution code
void exec();
+ void initizeCutOffs(const double RCut, const double WCut);
void examineInput(API::MatrixWorkspace_const_sptr binAdj, API::MatrixWorkspace_const_sptr detectAdj);
API::MatrixWorkspace_sptr setUpOutputWorkspace(const std::vector<double> & binParams, const API::SpectraDetectorMap * const specMap) const;
//these are the steps that are run on each individual spectrum
- void convertWavetoQ(const size_t specIndex, const bool doGravity, MantidVec & Qs) const;
+ size_t waveLengthCutOff(const size_t specInd) const;
+ void calculateNormalization(const size_t wavStart, const size_t specInd, API::MatrixWorkspace_const_sptr pixelAdj, double const * const binNorms, double const * const binNormEs, const MantidVec::iterator norm, const MantidVec::iterator normETo2) const;
void pixelWeight(API::MatrixWorkspace_const_sptr pixelAdj, const size_t specIndex, double & weight, double & error) const;
- void addWaveAdj(API::MatrixWorkspace_const_sptr pixelAdj, const MantidVec * const binNorms, const MantidVec * const binNormEs, MantidVec & outNorms, MantidVec & outETo2) const;
- void normToBinWidth(const size_t specIndex, const MantidVec & QIns, MantidVec & theNorms, MantidVec & errorSquared) const;
+ void addWaveAdj(const double * c, const double * Dc, MantidVec::iterator bInOut, MantidVec::iterator e2InOut) const;
+ void normToBinWidth(const size_t offSet, const size_t specIndex, const MantidVec::iterator theNorms, const MantidVec::iterator errorSquared) const;
+ void convertWavetoQ(const size_t specInd, const bool doGravity, const size_t offset, MantidVec::iterator Qs) const;
+ void getQBinPlus1(const MantidVec & OutQs, const double QToFind, MantidVec::const_iterator & loc) const;
void updateSpecMap(const size_t specIndex, API::SpectraDetectorMap * const specMap, const Geometry::ISpectraDetectorMap & inSpecMap, API::MatrixWorkspace_sptr outputWS) const;
+ void normalize(const MantidVec & normSum, const MantidVec & normError2, MantidVec & YOut, MantidVec & errors) const;
};
} // namespace Algorithms
diff --git a/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp b/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
index 7d242c9f424..329c8112fc1 100644
--- a/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
@@ -17,7 +17,7 @@ namespace Algorithms
{
// Register the algorithm into the AlgorithmFactory
-//DECLARE_ALGORITHM(Q1D2)
+DECLARE_ALGORITHM(Q1D2)
/// Sets documentation strings for this algorithm
void Q1D2::initDocs()
@@ -39,14 +39,21 @@ void Q1D2::init()
dataVal->add(new InstrumentValidator<>);
dataVal->add(new CommonBinsValidator<>);
declareProperty(new WorkspaceProperty<>("DetBankWorkspace", "", Direction::Input, dataVal));
+ declareProperty(new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output));
+ declareProperty(new ArrayProperty<double>("OutputBinning", new RebinParamsValidator));
+
declareProperty(new WorkspaceProperty<>("PixelAdj","", Direction::Input, true));
CompositeValidator<> *wavVal = new CompositeValidator<>;
wavVal->add(new WorkspaceUnitValidator<>("Wavelength"));
wavVal->add(new HistogramValidator<>);
declareProperty(new WorkspaceProperty<>("WavelengthAdj", "", Direction::Input, true, wavVal));
- declareProperty(new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output));
- declareProperty(new ArrayProperty<double>("OutputBinning", new RebinParamsValidator));
+
declareProperty("AccountForGravity",false);
+
+ BoundedValidator<double> *mustBePositive = new BoundedValidator<double>();
+ mustBePositive->setLower(0.0);
+ declareProperty("RadiusCut", 0.0, mustBePositive);
+ declareProperty("WaveCut", 0.0, mustBePositive->clone());
}
/**
@ throw invalid_argument if the workspaces are not mututially compatible
@@ -58,13 +65,14 @@ void Q1D2::exec()
// this pointer could be NULL as PixelAdj is an optional property
MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
const bool doGravity = getProperty("AccountForGravity");
+ initizeCutOffs(getProperty("RadiusCut"), getProperty("WaveCut"));
//throws if we don't have common binning or another incompatibility
examineInput(waveAdj, pixelAdj);
// normalization as a function of wavelength (i.e. centers of x-value bins)
- MantidVec const * const binNorms = waveAdj ? &(waveAdj->readY(0)) : NULL;
+ double const * const binNorms = waveAdj ? &(waveAdj->readY(0)[0]) : NULL;
// error on the wavelength normalization
- MantidVec const * const binNormEs = waveAdj ? &(waveAdj->readE(0)) : NULL;
+ double const * const binNormEs = waveAdj ? &(waveAdj->readE(0)[0]) : NULL;
//define the (large number of) data objects that are going to be used in all iterations of the loop below
// Construct a new spectra map. This will be faster than remapping the old one
@@ -80,9 +88,6 @@ void Q1D2::exec()
MantidVec normError2(YOut.size(), 0.0);
const Geometry::ISpectraDetectorMap & inSpecMap = m_dataWS->spectraMap();
- //const Axis* const spectraAxis = m_dataWS->getAxis(1);
-
-
const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
Progress progress(this, 0.1, 1.0, numSpec+1);
@@ -107,42 +112,57 @@ void Q1D2::exec()
continue;
}
- // A temporary vector to store the Q values for the input workspace before the rebin
- MantidVec QIn(m_numInBins);
- convertWavetoQ(i, doGravity, QIn);
- const MantidVec & YIn = m_dataWS->readY(i);
- const MantidVec & EIn = m_dataWS->readE(i);
-
- double detectorAdj, detAdjErr;
- pixelWeight(pixelAdj, i, detectorAdj, detAdjErr);
+ //get the bins that are included inside the RadiusCut/WaveCutcut off, those to calculate for
+ size_t wavStart(0);
+ if (m_RCut > 1e-200)
+ {
+ wavStart = waveLengthCutOff(i);
+ if (wavStart >= m_dataWS->readY(i).size())
+ {
+ // all the spectra in this detector is out of range
+ continue;
+ }
+ }
+
+ const size_t numWavbins = m_dataWS->readY(i).size()-wavStart;
+ //make just one call to new to reduce CPU overhead on each thread, access is these three arrays is via iterators
+ MantidVec _noDirectUseStorage_(3*numWavbins);
+ //normalization term
+ MantidVec::iterator norms = _noDirectUseStorage_.begin();
+ // the error on these weights, it contributes to the error calculation on the output workspace
+ MantidVec::iterator normETo2s = norms + numWavbins;
+ // the Q values calculated from input wavelength workspace
+ MantidVec::iterator QIn = normETo2s + numWavbins;
// the weighting for this input spectrum that is added to the normalization
- MantidVec norm(m_numInBins, detectorAdj);
- // the error on these weights, it contributes to the error calculation on the output workspace
- MantidVec normETo2(m_numInBins, detAdjErr*detAdjErr);
- addWaveAdj(pixelAdj, binNorms, binNormEs, norm, normETo2);
- normToBinWidth(i, QIn, norm, normETo2);
-
- //find the output bin that each input y-value will fall into, remembering there is one more bin boundary than bins
+ calculateNormalization(wavStart, i, pixelAdj, binNorms, binNormEs, norms, normETo2s);
+
+ // now read the data from the input workspace, calculate Q for each bin
+ convertWavetoQ(i, doGravity, wavStart, QIn);
+ // Pointers to the counts data and it's error
+ MantidVec::const_iterator YIn = m_dataWS->readY(i).begin()+wavStart;
+ MantidVec::const_iterator EIn = m_dataWS->readE(i).begin()+wavStart;
+
+ //when finding the output Q bin remember that the input Q bins (from the convert to wavelength) start high and reduce
MantidVec::const_iterator loc = QOut.end();
- for(size_t j=0; j < m_numInBins; ++j)
+ // sum the Q contributions from each individual spectrum into the output array
+ const MantidVec::const_iterator end = m_dataWS->readY(i).end();
+ for( ; YIn != end; ++YIn, ++EIn, ++QIn, ++norms, ++normETo2s)
{
- //Q goes from a high value to a low one in the QIn array (high Q particles arrive at low TOF) so we know loc will go downwards
- loc = std::upper_bound(QOut.begin(), loc, QIn[j]);
+ //find the output bin that each input y-value will fall into, remembering there is one more bin boundary than bins
+ getQBinPlus1(QOut, *QIn, loc);
// ignore counts that are out of the output range
if ( (loc != QOut.begin()) && (loc != QOut.end()) )
{
const size_t bin = loc - QOut.begin() - 1;
PARALLEL_CRITICAL(q1d_counts_sum)
{
- YOut[bin] += YIn[j];
- normSum[bin] += norm[j];
+ YOut[bin] += *YIn;
+ normSum[bin] += *norms;
//these are the errors squared which will be summed and square rooted at the end
- EOutTo2[bin] += EIn[j]*EIn[j];
- normError2[bin] += normETo2[j];
+ EOutTo2[bin] += (*EIn)*(*EIn);
+ normError2[bin] += *normETo2s;
}
- //this is used to restrict the search range above for a modest increase in speed
- ++loc;
}
}
@@ -156,19 +176,9 @@ void Q1D2::exec()
}
PARALLEL_CHECK_INTERUPT_REGION
- for (size_t k = 0; k < YOut.size(); ++k)
- {
- // the normalisation is a = b/c where b = counts c =normalistion term
- const double c = normSum[k];
- const double a = YOut[k] /= c;
- // when a = b/c, the formula for Da, the error on a, in terms of Db, etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
- //(Da)^2 = ((Db/b)^2 + (Dc/c)^2)*(b^2/c^2) = ((Db/c)^2 + (b*Dc/c^2)^2) = (Db^2 + (b*Dc/c)^2)/c^2 = (Db^2 + (Dc*a)^2)/c^2
- //this will work as long as c>0, but then the above formula above can't deal with 0 either
- const double aOverc = a/c;
- EOutTo2[k] = std::sqrt(EOutTo2[k]/(c*c) + normError2[k]*aOverc*aOverc);
-
- progress.report("Computing I(Q)");
- }
+ progress.report("Normalizing I(Q)");
+ //finally divide the number of counts in each output Q bin by its weighting
+ normalize(normSum, normError2, YOut, EOutTo2);
setProperty("OutputWorkspace",outputWS);
}
@@ -184,7 +194,6 @@ void Q1D2::examineInput(API::MatrixWorkspace_const_sptr binAdj, API::MatrixWorks
{
throw std::invalid_argument("Empty data workspace passed, can not continue");
}
- m_numInBins = m_dataWS->readY(0).size();
//it is not an error for these workspaces not to exist
if (binAdj)
@@ -193,7 +202,7 @@ void Q1D2::examineInput(API::MatrixWorkspace_const_sptr binAdj, API::MatrixWorks
{
throw std::invalid_argument("The WavelengthAdj workspace must have one spectrum");
}
- if ( binAdj->readY(0).size() != m_numInBins )
+ if ( binAdj->readY(0).size() != m_dataWS->readY(0).size() )
{
throw std::invalid_argument("The WavelengthAdj workspace's bins must match those of the detector bank workspace");
}
@@ -232,6 +241,18 @@ if ( binAdj->isDistribution() != m_dataWS->isDistribution() )
g_log.debug() << "All input workspaces were found to be valid\n";
}
+/** Detector independent parts of the wavelength cut off calculation
+* @param RCut the radius cut off, should be value of the property RadiusCut
+* @param WCut this wavelength cut off, should be equal to the value WaveCut
+*/
+void Q1D2::initizeCutOffs(const double RCut, const double WCut)
+{
+ if ( RCut > 0 && WCut > 0 )
+ {
+ m_WCutOver = WCut/RCut;
+ m_RCut = RCut;
+ }
+}
/** Creates the output workspace, its size, units, etc.
* @param binParams the bin boundary specification using the same same syntax as param the Rebin algorithm
* @param specMap a spectra map that the new workspace should use and take owner ship of
@@ -256,46 +277,52 @@ API::MatrixWorkspace_sptr Q1D2::setUpOutputWorkspace(const std::vector<double> &
outputWS->replaceSpectraMap(specMap);
return outputWS;
}
-/** Fills a vector with the Q values calculated from the wavelengths in the input workspace and the workspace
-* geometry as Q = 4*pi*sin(theta)/lambda
-* @param[in] specIndex the spectrum to calculate
-* @param[in] doGravity if to include gravity in the calculation of Q
-* @param[out] Qs a preallocated array that is large enough to contain all the calculated Q values
-* @throw NotFoundError if the detector associated with the spectrum is not found in the instrument definition
+/** Finds the first index number of the first wavelength bin that should included based on the
+* the calculation: W = Wcut (Rcut-R)/Rcut. Must only be called if m_RCut > 0.0
+* @param specInd spectrum that is being analysed
+* @return index number of the first bin to include in the calculation
*/
-void Q1D2::convertWavetoQ(const size_t specIndex, const bool doGravity, MantidVec & Qs) const
+size_t Q1D2::waveLengthCutOff(const size_t specInd) const
{
- static const double FOUR_PI=4.0*M_PI;
- const MantidVec & XIn = m_dataWS->readX(specIndex);
- IDetector_const_sptr det = m_dataWS->getDetector(specIndex);
+ //get the distance of between this detector and the origin, which should be the along the beam center
+ const V3D posOnBank = m_dataWS->getDetector(specInd)->getPos();
+ double R = (posOnBank.X()*posOnBank.X())+(posOnBank.Y()*posOnBank.Y());
+ R = std::sqrt(R);
- if (doGravity)
+ const double WMin = m_WCutOver*(m_RCut-R);
+ const MantidVec & Xs = m_dataWS->readX(specInd);
+ return std::lower_bound(Xs.begin(), Xs.end(), WMin) - Xs.begin();
+}
+/** Calcualtes the normalization term for each output bin
+* @param[in] offSet the inex number of the first bin in the input wavelengths that is actually being used
+* @param[in] specInd the spectrum to calculate
+* @param[in] pixelAdj if not NULL this is workspace contains single bins with the adjustments, e.g. detector efficencies, for the given spectrum index
+* @param[in] binNorms pointer to a contigious array of doubles that are the wavelength correction from waveAdj workspace, can be NULL
+* @param[in] binNormEs pointer to a contigious array of doubles which corrospond to the corrections and are their errors, can be NULL
+* @param[out] norm normalization for each bin, including soild angle, pixel correction, the proportion that is not masked and the normalization workspace
+* @param[out] normETo2 this pointer must point to the end of the norm array, it will be filled with the total of the error on the normalization
+*/
+void Q1D2::calculateNormalization(const size_t wavStart, const size_t specInd, API::MatrixWorkspace_const_sptr pixelAdj, double const * const binNorms, double const * const binNormEs, const MantidVec::iterator norm, const MantidVec::iterator normETo2) const
+{
+ double detectorAdj, detAdjErr;
+ pixelWeight(pixelAdj, specInd, detectorAdj, detAdjErr);
+
+ //use that the normalization array ends at the start of the error array
+ for(MantidVec::iterator n = norm, e = normETo2; n != normETo2; ++n, ++e)
{
- GravitySANSHelper grav(m_dataWS, det);
- for ( size_t j = 0; j < m_numInBins; ++j)
- {
- // as the fall under gravity is wavelength dependent sin theta is now different for each bin with each detector
- // the HistogramValidator at the start should ensure that we have one more bin on the readX()s
- const double lambda = (XIn[j]+XIn[j+1])/2.0;
- const double sinTheta = grav.calcSinTheta(lambda);
- // Now we're ready to go to Q
- Qs[j] = FOUR_PI*sinTheta/lambda;
- }
+ *n = detectorAdj;
+ *e = detAdjErr*detAdjErr;
}
- else
+
+ if (binNorms && binNormEs)
{
- // Calculate the Q values for the current spectrum, using Q = 4*pi*sin(theta)/lambda
- const double factor = 2.0* FOUR_PI*sin( m_dataWS->detectorTwoTheta(det)/2.0 );
- for ( size_t j = 0; j < m_numInBins; ++j)
- {
- // the HistogramValidator at the start should ensure that we have one more bin on the readX()s
- Qs[j] = factor/(XIn[j]+XIn[j+1]);
- }
+ addWaveAdj(binNorms+wavStart, binNormEs+wavStart, norm, normETo2);
}
+ normToBinWidth(wavStart, specInd, norm, normETo2);
}
/** Calculates the normalisation for the spectrum specified by the index number that was passed
* as the solid anlge multiplied by the pixelAdj that was passed
-* @param[in] pixelAdj if not NULL this is workspace contains single bins with the adjustments, e.g. detector efficencie, for hte given spectrum index
+* @param[in] pixelAdj if not NULL this is workspace contains single bins with the adjustments, e.g. detector efficencies, for the given spectrum index
* @param[in] specIndex the spectrum index to return the data from
* @param[out] weight the solid angle or if pixelAdj the solid anlge times the pixel adjustment for this spectrum
* @param[out] error the error on the weight, only non-zero if pixelAdj
@@ -322,50 +349,44 @@ void Q1D2::pixelWeight(API::MatrixWorkspace_const_sptr pixelAdj, const size_t s
}
}
/** Calculates the contribution to the normalization terms from each bin in a spectrum
-* @param[in] pixelAdj detector efficiency input workspace
-* @param[in] binNorms the wavelength dependent normalization term
-* @param[in] binNormEs the wavelength dependent term's error
-* @param[in,out] outNorms normalization for each bin, this method multiplise this by the proportion that is not masked and the normalization workspace
-* @param[in, out] outETo2 the error on the normalisation term before the WavelengthAdj term
+* @param[in] c pointer to the start of a contigious array of wavelength dependent normalization terms
+* @param[in] Dc pointer to the start of a contigious array that corrosponds to wavelength dependent term, having its error
+* @param[in,out] bInOut normalization for each bin, this method multiplise this by the proportion that is not masked and the normalization workspace
+* @param[in, out] e2InOut this array must follow straight after the normalization array and will contain the error on the normalisation term before the WavelengthAdj term
*/
-void Q1D2::addWaveAdj(API::MatrixWorkspace_const_sptr pixelAdj, const MantidVec * const binNorms, const MantidVec * const binNormEs, MantidVec & outNorms, MantidVec & outETo2) const
+void Q1D2::addWaveAdj(const double * c, const double * Dc, MantidVec::iterator bInOut, MantidVec::iterator e2InOut) const
{
- UNUSED_ARG(pixelAdj)
// normalize by the wavelength dependent correction, keeping the percentage errors the same
- if (binNorms && binNormEs)
+ // the error when a = b*c, the formula for Da, the error on a, in terms of Db, etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
+ //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) = (Db*c)^2 + (Dc*b)^2
+ // the variable names relate to those above as: existing values (b=bInOut) multiplied by the additional errors (Dc=binNormEs), existing errors (Db=sqrt(e2InOut)) times new factor (c=binNorms)
+
+ //use the fact that errors array follows straight after the normalization array
+ const MantidVec::const_iterator end = e2InOut;
+ for( ; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut)
{
- // the error when a = b*c, the formula for Da, the error on a, in terms of Db, etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
- //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) = (Db*c)^2 + (Dc*b)^2
- // name the things to fit with the equation above: existing values (b=bInOut(=outNorms)) multiplied by the additional errors (Dc=binNormEs), existing errors (Db=sqrt(e2InOut(=outETo2))) times new factor (c=binNorms)
- MantidVec::iterator e2InOut=outETo2.begin(), bInOut=outNorms.begin();
- MantidVec::const_iterator c=binNorms->begin(), Dc=binNormEs->begin();
- MantidVec::const_iterator end = outETo2.end();
- for( ; e2InOut != end; ++e2InOut, ++c, ++Dc, ++bInOut)
- {
- //first the error
- *e2InOut += ( (*e2InOut)*(*c)*(*c) )+( (*Dc)*(*Dc)*(*bInOut)*(*bInOut) );
- // now the actual calculation a = b*c
- *bInOut = (*bInOut)*(*c);
- }
+ //first the error
+ *e2InOut += ( (*e2InOut)*(*c)*(*c) )+( (*Dc)*(*Dc)*(*bInOut)*(*bInOut) );
+ // now the actual calculation a = b*c
+ *bInOut = (*bInOut)*(*c);
}
}
/** Add the bin widths, scaled to bin masking, to the normalization
+* @param[in] offSet the inex number of the first bin in the input wavelengths that is actually being used
* @param[in] specIndex the spectrum to calculate
-* @param[in] QIns Q bin boundaries
* @param[in,out] theNorms normalization for each bin, this is multiplied by the proportion that is not masked and the normalization workspace
-* @param[in,out] errorSquared the running total of the error on the normalization
+* @param[in,out] errorSquared the running total of the square of the uncertainty in the normalization
*/
-void Q1D2::normToBinWidth(const size_t specIndex, const MantidVec & QIns, MantidVec & theNorms, MantidVec & errorSquared) const
+void Q1D2::normToBinWidth(const size_t offSet, const size_t specIndex, const MantidVec::iterator theNorms, const MantidVec::iterator errorSquared) const
{
- UNUSED_ARG(QIns)
/* //normally this is false but handling this would mean more combinations of distribution/raw counts workspaces could be accepted
if (m_convToDistr)
{
for(int i = 0; i < theNorms.size(); ++i)
{
const double width = ???;
- theNorms[i] *= width;
- errorSquared[i] *= width*width;
+ *(theNorms+i) *= width;
+ *(errorSquared+i) *= width*width;
}
}*/
@@ -378,14 +399,99 @@ void Q1D2::normToBinWidth(const size_t specIndex, const MantidVec & QIns, Mantid
MatrixWorkspace::MaskList::const_iterator it;
for (it = mask.begin(); it != mask.end(); ++it)
{
+ size_t outBin = it->first;
+ if ( outBin < offSet )
+ {
+ // this masked bin wasn't in the range being delt with anyway
+ continue;
+ }
+ outBin -= offSet;
// The weight for this masked bin is 1 - the degree to which this bin is masked
const double factor = 1.0-(it->second);
- theNorms[it->first] *= factor;
- errorSquared[it->first] *= factor*factor;
+ *(theNorms+outBin) *= factor;
+ *(errorSquared+outBin) *= factor*factor;
}
}
}
-/** !!!PROTOTYPE needs more testing !!! Map all the detectors onto the spectrum of the output
+/** Fills a vector with the Q values calculated from the wavelength bin centers from the input workspace and
+* the workspace geometry as Q = 4*pi*sin(theta)/lambda
+* @param[in] specInd the spectrum to calculate
+* @param[in] doGravity if to include gravity in the calculation of Q
+* @param[in] offset index number of the first input bin to use
+* @param[out] Qs points to a preallocated array that is large enough to contain all the calculated Q values
+* @throw NotFoundError if the detector associated with the spectrum is not found in the instrument definition
+*/
+void Q1D2::convertWavetoQ(const size_t specInd, const bool doGravity, const size_t offset, MantidVec::iterator Qs) const
+{
+ static const double FOUR_PI=4.0*M_PI;
+
+ IDetector_const_sptr det = m_dataWS->getDetector(specInd);
+
+ // wavelengths (lamda) to be converted to Q
+ MantidVec::const_iterator waves = m_dataWS->readX(specInd).begin() + offset;
+ // going from bin boundaries to bin centered x-values the size goes down one
+ const MantidVec::const_iterator end = m_dataWS->readX(specInd).end() - 1;
+ if (doGravity)
+ {
+ GravitySANSHelper grav(m_dataWS, det);
+ for( ; waves !=end; ++Qs, ++waves)
+ {
+ // the HistogramValidator at the start should ensure that we have one more bin on the input wavelengths
+ const double lambda = 0.5*(*(waves+1)+(*waves));
+ // as the fall under gravity is wavelength dependent sin theta is now different for each bin with each detector
+ const double sinTheta = grav.calcSinTheta(lambda);
+ // Now we're ready to go to Q
+ *Qs = FOUR_PI*sinTheta/lambda;
+ }
+ }
+ else
+ {
+ // Calculate the Q values for the current spectrum, using Q = 4*pi*sin(theta)/lambda
+ const double factor = 2.0* FOUR_PI*sin( m_dataWS->detectorTwoTheta(det)/2.0 );
+ for( ; waves !=end; ++Qs, ++waves)
+ {
+ // the HistogramValidator at the start should ensure that we have one more bin on the input wavelengths
+ *Qs = factor/(*(waves+1)+(*waves));
+ }
+ }
+}
+/** This is a slightly "clever" method as it makes some guesses about where is best
+* to look for the right Q bin based on the fact that the input Qs (calcualted from wavelengths) tend
+* to go down while the output Qs are always in accending order
+* @param[in] OutQs the array of output Q bin boundaries, this finds the bin that contains the QIn value
+* @param[in] QToFind the Q value to find the correct bin for
+* @param[in, out] loc points to the bin boundary (in the OutQs array) whos Q is higher than QToFind and higher by the smallest amount. Algorithm starts by checking the value of loc passed and then all the bins _downwards_ through the array
+*/
+void Q1D2::getQBinPlus1(const MantidVec & OutQs, const double QToFind, MantidVec::const_iterator & loc) const
+{
+ if ( loc != OutQs.end() )
+ {
+ while ( loc != OutQs.begin() )
+ {
+ if ( (QToFind >= *(loc-1)) && (QToFind < *loc) )
+ {
+ return;
+ }
+ --loc;
+ }
+ if ( QToFind < *loc )
+ {
+ //QToFind was outside the array
+ return;
+ }
+ }
+ else //loc == OutQs.end()
+ {
+ if ( OutQs.empty() || QToFind > *(loc-1) )
+ {
+ return;
+ }
+ }
+
+ // we are lost, normally the order of the Q values means we only get here on the first iteration. It's slow
+ loc = std::lower_bound(OutQs.begin(), OutQs.end(), QToFind);
+}
+/** Map all the detectors onto the spectrum of the output
* @param[in] specIndex the spectrum to add
* @param[out] specMap the map in the output workspace to write to
* @param[in] inSpecMap spectrum data
@@ -400,6 +506,27 @@ void Q1D2::updateSpecMap(const size_t specIndex, API::SpectraDetectorMap * const
specMap->addSpectrumEntries(newSpectrumNo,inSpecMap.getDetectors(spectraAxis->spectraNo(specIndex)));
}
}
+/** Divides the number of counts in each output Q bin by the wrighting ("number that would expected to arrive")
+* The errors are propogated using the uncorrolated error estimate for multiplication/division
+* @param[in] normSum the weighting for each bin
+* @param[in] normError2 square of the error on the normalization
+* @param[in, out] counts counts in each bin
+* @param[in, out] errors input the _square_ of the error on each bin, output the total error (unsquared)
+*/
+void Q1D2::normalize(const MantidVec & normSum, const MantidVec & normError2, MantidVec & counts, MantidVec & errors) const
+{
+ for (size_t k = 0; k < counts.size(); ++k)
+ {
+ // the normalisation is a = b/c where b = counts c =normalistion term
+ const double c = normSum[k];
+ const double a = counts[k] /= c;
+ // when a = b/c, the formula for Da, the error on a, in terms of Db, etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
+ //(Da)^2 = ((Db/b)^2 + (Dc/c)^2)*(b^2/c^2) = ((Db/c)^2 + (b*Dc/c^2)^2) = (Db^2 + (b*Dc/c)^2)/c^2 = (Db^2 + (Dc*a)^2)/c^2
+ //this will work as long as c>0, but then the above formula above can't deal with 0 either
+ const double aOverc = a/c;
+ errors[k] = std::sqrt(errors[k]/(c*c) + normError2[k]*aOverc*aOverc);
+ }
+}
} // namespace Algorithms
} // namespace Mantid
diff --git a/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h b/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
index 63d3640290b..80b3782c12e 100644
--- a/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
+++ b/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
@@ -114,7 +114,7 @@ public:
Mantid::API::AnalysisDataService::Instance().remove(m_noGrav);
}
-
+
void testGravity()
{
Mantid::Algorithms::Q1D2 Q1D;
@@ -157,6 +157,48 @@ public:
Mantid::API::AnalysisDataService::Instance().remove(outputWS);
}
+
+ void testCuts()
+ {
+ Mantid::Algorithms::Q1D2 Q1D;
+ Q1D.initialize();
+
+ TS_ASSERT_THROWS_NOTHING(
+ Q1D.setProperty("DetBankWorkspace", m_inputWS);
+ Q1D.setProperty("WavelengthAdj", m_wavNorm);
+ Q1D.setPropertyValue("PixelAdj", m_pixel);
+ Q1D.setPropertyValue("OutputWorkspace", m_noGrav);
+ Q1D.setPropertyValue("OutputBinning", "0.1,-0.02,0.5");
+ Q1D.setProperty("RadiusCut", 0.22);
+ Q1D.setProperty("WaveCut", 8.0);
+ )
+ //#default is don't correct for gravity
+ TS_ASSERT_THROWS_NOTHING( Q1D.execute() )
+ TS_ASSERT( Q1D.isExecuted() )
+
+ Mantid::API::MatrixWorkspace_sptr result;
+ TS_ASSERT_THROWS_NOTHING( result = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::AnalysisDataService::Instance().retrieve(m_noGrav)) )
+ TS_ASSERT(result)
+ TS_ASSERT_EQUALS( result->getNumberHistograms(), 1 )
+
+ TS_ASSERT_EQUALS( result->readX(0).size(), 83 )
+ TS_ASSERT_EQUALS( result->readX(0).front(), 0.1 )
+ TS_ASSERT_DELTA( result->readX(0)[3], 0.1061208, 1e-6 )
+ TS_ASSERT_DELTA( result->readX(0)[56], 0.3031165, 1e-5 )
+ TS_ASSERT_EQUALS( result->readX(0).back(), 0.5 )
+
+ TS_ASSERT_DELTA( result->readY(0).front(), 1192471.95, 0.1 )
+ TS_ASSERT( boost::math::isnan(result->readY(0)[3]))
+ TS_ASSERT_DELTA( result->readY(0)[12], 503242.79, 0.1)
+ TS_ASSERT( boost::math::isnan(result->readY(0).back()) )
+
+ TS_ASSERT_DELTA( result->readE(0)[2], 4847257060, 10 )
+ TS_ASSERT_DELTA( result->readE(0)[10], 4921866100, 100 )
+ TS_ASSERT( boost::math::isnan(result->readE(0)[7]) )
+
+ Mantid::API::AnalysisDataService::Instance().remove(m_noGrav);
+ }
void testInvalidInput()
{
--
GitLab