From b04a8108ea15e21f97d70aa98875647fb63d5890 Mon Sep 17 00:00:00 2001
From: Alex Buts <Alex.Buts@stfc.ac.uk>
Date: Tue, 27 Nov 2012 15:29:12 +0000
Subject: [PATCH] refs #6230 Modified sum spectra to produce weighted sum on
request
---
.gitignore | 1 +
.../inc/MantidAlgorithms/SumSpectra.h | 9 +-
.../Framework/Algorithms/src/SumSpectra.cpp | 130 ++++++++++++++----
.../Algorithms/test/SumSpectraTest.h | 114 ++++++++++++++-
4 files changed, 221 insertions(+), 33 deletions(-)
diff --git a/.gitignore b/.gitignore
index 3a6cdc58602..dfee419542c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -163,3 +163,4 @@ Desktop.ini
# Mac OS X Finder
.DS_Store
+/Code\Mantid\Framework\API\test/*.cpp
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SumSpectra.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SumSpectra.h
index ecf4a074a83..a0064014fac 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SumSpectra.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SumSpectra.h
@@ -69,10 +69,12 @@ public:
private:
/// Handle logic for RebinnedOutput workspaces
void doRebinnedOutput(API::MatrixWorkspace_sptr outputWorkspace,
- API::Progress &progress);
+ API::Progress &progress,
+ size_t &numSpectra,size_t &numMasked,size_t &numZeros);
/// Handle logic for Workspace2D workspaces
void doWorkspace2D(API::MatrixWorkspace_const_sptr localworkspace,
- API::ISpectrum *outSpec, API::Progress &progress);
+ API::ISpectrum *outSpec, API::Progress &progress,
+ size_t &numSpectra,size_t &numMasked,size_t &numZeros);
/// Sets documentation strings for this algorithm
virtual void initDocs();
// Overridden Algorithm methods
@@ -95,6 +97,9 @@ private:
int yLength;
/// Set of indicies to sum
std::set<int> indices;
+
+ // if calculateing additional workspace with specially weighted averages is necessary
+ bool m_CalculateWeightedSum;
};
} // namespace Algorithm
diff --git a/Code/Mantid/Framework/Algorithms/src/SumSpectra.cpp b/Code/Mantid/Framework/Algorithms/src/SumSpectra.cpp
index 90b728eba19..1dc5cbdda2a 100644
--- a/Code/Mantid/Framework/Algorithms/src/SumSpectra.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/SumSpectra.cpp
@@ -60,6 +60,10 @@ void SumSpectra::init()
"Can be specified instead of in addition to StartWorkspaceIndex and EndWorkspaceIndex.");
declareProperty("IncludeMonitors",true,"Whether to include monitor spectra in the sum (default: yes)");
+
+ declareProperty("WeightedSum",false,"Instead of the usual spectra sum, calculate the weighted sum in the form: \n"
+ "<math>\\Sigma(Signal_i/Error_i^2)/\\Sigma(1/Error_i^2)</math> This property is ignored for event workspace.\n"
+ "The sums are defined for <math>Error_i != 0</math> only, so the values are NaN if error for some value is equal to 0.");
}
/** Executes the algorithm
@@ -115,9 +119,13 @@ void SumSpectra::exec()
g_log.information() << "Spectra remapping gives single spectra with spectra number: "
<< m_outSpecId << "\n";
+ m_CalculateWeightedSum = getProperty("WeightedSum");
+
+
EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast<const EventWorkspace>(localworkspace);
if (eventW)
{
+ m_CalculateWeightedSum = false;
this->execEvent(eventW, this->indices);
}
else
@@ -127,7 +135,10 @@ void SumSpectra::exec()
// Create the 2D workspace for the output
MatrixWorkspace_sptr outputWorkspace = API::WorkspaceFactory::Instance().create(localworkspace,
1,localworkspace->readX(0).size(),this->yLength);
-
+ size_t numSpectra(0); // total number of processed spectra
+ size_t numMasked(0); // total number of the masked and skipped spectra
+ size_t numZeros(0); // number of spectra which have 0 value in the first column (used in special cases of evaluating how good Puasonian statistics is)
+
Progress progress(this, 0, 1, this->indices.size());
// This is the (only) output spectrum
@@ -142,11 +153,11 @@ void SumSpectra::exec()
if (localworkspace->id() == "RebinnedOutput")
{
- this->doRebinnedOutput(outputWorkspace, progress);
+ this->doRebinnedOutput(outputWorkspace, progress,numSpectra,numMasked,numZeros);
}
else
{
- this->doWorkspace2D(localworkspace, outSpec, progress);
+ this->doWorkspace2D(localworkspace, outSpec, progress,numSpectra,numMasked,numZeros);
}
// Pointer to sqrt function
@@ -157,6 +168,11 @@ void SumSpectra::exec()
std::transform(YError.begin(), YError.end(), YError.begin(), rs);
outputWorkspace->generateSpectraMap();
+ // set up the summing statistics
+ outputWorkspace->mutableRun().addProperty("NumAllSpectra",numSpectra,"",true);
+ outputWorkspace->mutableRun().addProperty("NumMaskSpectra",numMasked,"",true);
+ outputWorkspace->mutableRun().addProperty("NumZeroSpectra",numZeros,"",true);
+
// Assign it to the output workspace property
setProperty("OutputWorkspace", outputWorkspace);
@@ -199,12 +215,17 @@ specid_t SumSpectra::getOutputSpecId(MatrixWorkspace_const_sptr localworkspace)
* @param progress the progress indicator
*/
void SumSpectra::doWorkspace2D(MatrixWorkspace_const_sptr localworkspace,
- ISpectrum *outSpec, Progress &progress)
+ ISpectrum *outSpec, Progress &progress,
+ size_t &numSpectra,size_t &numMasked,size_t &numZeros)
{
// Get references to the output workspaces's data vectors
MantidVec& YSum = outSpec->dataY();
MantidVec& YError = outSpec->dataE();
+ MantidVec Weight;
+ if(m_CalculateWeightedSum)Weight.assign(YSum.size(),0);
+
+
// Loop over spectra
std::set<int>::iterator it;
//for (int i = m_MinSpec; i <= m_MaxSpec; ++i)
@@ -225,28 +246,60 @@ void SumSpectra::doWorkspace2D(MatrixWorkspace_const_sptr localworkspace,
// Skip monitors, if the property is set to do so
if ( !keepMonitors && det->isMonitor() ) continue;
// Skip masked detectors
- if ( det->isMasked() ) continue;
+ if ( det->isMasked() )
+ {
+ numMasked++;
+ continue;
+ }
}
catch(...)
{
// if the detector not found just carry on
}
+ numSpectra++;
// Retrieve the spectrum into a vector
const MantidVec& YValues = localworkspace->readY(i);
const MantidVec& YErrors = localworkspace->readE(i);
-
- for (int k = 0; k < this->yLength; ++k)
+ if(m_CalculateWeightedSum)
{
- YSum[k] += YValues[k];
- YError[k] += YErrors[k]*YErrors[k];
+ for (int k = 0; k < this->yLength; ++k)
+ {
+ if(YErrors[k]!=0)
+ {
+ double errsq = YErrors[k]*YErrors[k];
+ YError[k] +=errsq;
+ Weight[k] +=1./errsq;
+ YSum[k] += YValues[k]/errsq;
+ }
+ else
+ {
+ YSum[k] =std::numeric_limits<double>::quiet_NaN();
+ Weight[k]=std::numeric_limits<double>::quiet_NaN();
+ }
+ }
+ }
+ else
+ {
+ for (int k = 0; k < this->yLength; ++k)
+ {
+ YSum[k] += YValues[k];
+ YError[k] += YErrors[k]*YErrors[k];
+ }
}
// Map all the detectors onto the spectrum of the output
outSpec->addDetectorIDs( localworkspace->getSpectrum(i)->getDetectorIDs() );
+ if(YSum[0]==0)numZeros++;
progress.report();
}
+ if(m_CalculateWeightedSum)
+ {
+ for(size_t i=0;i<Weight.size();i++)
+ YSum[i]/=Weight[i];
+ }
+
}
/**
@@ -255,7 +308,8 @@ void SumSpectra::doWorkspace2D(MatrixWorkspace_const_sptr localworkspace,
* @param progress the progress indicator
*/
void SumSpectra::doRebinnedOutput(MatrixWorkspace_sptr outputWorkspace,
- Progress &progress)
+ Progress &progress,
+ size_t &numSpectra,size_t &numMasked,size_t &numZeros)
{
// Get a copy of the input workspace
MatrixWorkspace_sptr temp = getProperty("InputWorkspace");
@@ -283,12 +337,11 @@ void SumSpectra::doRebinnedOutput(MatrixWorkspace_sptr outputWorkspace,
MantidVec& YSum = outSpec->dataY();
MantidVec& YError = outSpec->dataE();
MantidVec& FracSum = outWS->dataF(0);
+ MantidVec Weight;
+ if(m_CalculateWeightedSum)Weight.assign(YSum.size(),0);
// Loop over spectra
std::set<int>::iterator it;
- size_t numSpectra(0); // total number of processed spectra
- size_t numMasked(0); // total number of the masked and skipped spectra
- size_t numZeros(0); // number of spectra which have 0 value in the first column (used in special cases of evaluating how good Puasonian statistics is)
//for (int i = m_MinSpec; i <= m_MaxSpec; ++i)
for (it = indices.begin(); it != indices.end(); ++it)
{
@@ -300,7 +353,6 @@ void SumSpectra::doRebinnedOutput(MatrixWorkspace_sptr outputWorkspace,
break;
}
- numSpectra++;
try
{
// Get the detector object for this spectrum
@@ -318,19 +370,43 @@ void SumSpectra::doRebinnedOutput(MatrixWorkspace_sptr outputWorkspace,
{
// if the detector not found just carry on
}
+ numSpectra++;
// Retrieve the spectrum into a vector
const MantidVec& YValues = localworkspace->readY(i);
const MantidVec& YErrors = localworkspace->readE(i);
const MantidVec& FracArea = inWS->readF(i);
- for (int k = 0; k < this->yLength; ++k)
+ if(m_CalculateWeightedSum)
+ {
+ for (int k = 0; k < this->yLength; ++k)
+ {
+ if(YErrors[k]!=0)
+ {
+ double errsq = YErrors[k]*YErrors[k]*FracArea[k]*FracArea[k];
+ YError[k] +=errsq;
+ Weight[k] +=1./errsq;
+ YSum[k] += YValues[k]*FracArea[k]/errsq;
+ FracSum[k] += FracArea[k];
+ }
+ else
+ {
+ YSum[k] =std::numeric_limits<double>::quiet_NaN();
+ Weight[k]=std::numeric_limits<double>::quiet_NaN();
+ FracSum[k] += FracArea[k];
+ }
+ }
+ }
+ else
{
- YSum[k] += YValues[k] * FracArea[k];
- YError[k] += YErrors[k] * YErrors[k] * FracArea[k] * FracArea[k];
- FracSum[k] += FracArea[k];
+ for (int k = 0; k < this->yLength; ++k)
+ {
+ YSum[k] += YValues[k] * FracArea[k];
+ YError[k] += YErrors[k] * YErrors[k] * FracArea[k] * FracArea[k];
+ FracSum[k] += FracArea[k];
+ }
}
- if(YValues[0]==0)numZeros++;
+ if(YSum[0]==0)numZeros++;
// Map all the detectors onto the spectrum of the output
outSpec->addDetectorIDs(localworkspace->getSpectrum(i)->getDetectorIDs());
@@ -338,11 +414,7 @@ void SumSpectra::doRebinnedOutput(MatrixWorkspace_sptr outputWorkspace,
progress.report();
}
- outWS->mutableRun().addProperty("NumAllSpectra",numSpectra,"",true);
- outWS->mutableRun().addProperty("NumMaskSpectra",numMasked,"",true);
- outWS->mutableRun().addProperty("NumZeroSpectra",numZeros,"",true);
-
- // Create the correct representation
+ // Create the correct representation
outWS->finalize();
}
@@ -367,10 +439,9 @@ void SumSpectra::execEvent(EventWorkspace_const_sptr localworkspace, std::set<in
// Loop over spectra
std::set<int>::iterator it;
- size_t numSpectra(0); // total number of processed spectra
- size_t numMasked(0); // total number of the masked and skipped spectra
- size_t numZeros(0); // number of fully empty spectra with no events (used in special cases of evaluating how good Puasonian statistics is)
-
+ size_t numSpectra(0);
+ size_t numMasked(0);
+ size_t numZeros(0);
//for (int i = m_MinSpec; i <= m_MaxSpec; ++i)
for (it = indices.begin(); it != indices.end(); ++it)
{
@@ -382,7 +453,6 @@ void SumSpectra::execEvent(EventWorkspace_const_sptr localworkspace, std::set<in
break;
}
- numSpectra++;
try
{
// Get the detector object for this spectrum
@@ -400,6 +470,8 @@ void SumSpectra::execEvent(EventWorkspace_const_sptr localworkspace, std::set<in
{
// if the detector not found just carry on
}
+ numSpectra++;
+
//Add the event lists with the operator
const EventList & tOutEL = localworkspace->getEventList(i);
if(tOutEL.empty())
diff --git a/Code/Mantid/Framework/Algorithms/test/SumSpectraTest.h b/Code/Mantid/Framework/Algorithms/test/SumSpectraTest.h
index b53f7297c11..60ec923ce59 100644
--- a/Code/Mantid/Framework/Algorithms/test/SumSpectraTest.h
+++ b/Code/Mantid/Framework/Algorithms/test/SumSpectraTest.h
@@ -8,6 +8,7 @@
#include "MantidDataObjects/Workspace2D.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include <boost/lexical_cast.hpp>
#include <cxxtest/TestSuite.h>
using namespace Mantid;
@@ -22,8 +23,9 @@ public:
SumSpectraTest()
{
- inputSpace = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(10,102,true);
- inputSpace->instrumentParameters().addBool(inputSpace->getDetector(1).get(),"masked",true);
+ this->nTestHist=10;
+ this->inputSpace = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(nTestHist,102,true);
+ this->inputSpace->instrumentParameters().addBool(inputSpace->getDetector(1).get(),"masked",true);
}
~SumSpectraTest()
@@ -82,6 +84,15 @@ public:
TS_ASSERT_EQUALS( spec->getDetectorIDs().size(), 2);
TS_ASSERT( spec->hasDetectorID(3) );
TS_ASSERT( spec->hasDetectorID(4) );
+
+ TS_ASSERT(output2D->run().hasProperty("NumAllSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumMaskSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumZeroSpectra"))
+
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(3-1),output2D->run().getLogData("NumAllSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(1),output2D->run().getLogData("NumMaskSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(0),output2D->run().getLogData("NumZeroSpectra")->value())
+
}
void testExecWithoutLimits()
@@ -140,6 +151,17 @@ public:
TS_ASSERT( spec->hasDetectorID(6) );
TS_ASSERT( spec->hasDetectorID(7) );
+
+ TS_ASSERT(output2D->run().hasProperty("NumAllSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumMaskSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumZeroSpectra"))
+
+ // Spectra at workspace index 1 is masked, 8 & 9 are monitors
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(nTestHist-3),output2D->run().getLogData("NumAllSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(1),output2D->run().getLogData("NumMaskSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(0),output2D->run().getLogData("NumZeroSpectra")->value())
+
+
}
void testExecEvent_inplace()
@@ -186,6 +208,11 @@ public:
TS_ASSERT_EQUALS(output->getNumberHistograms(), 1);
TS_ASSERT_EQUALS(output->getNumberEvents(), 9 * numEvents);
TS_ASSERT_EQUALS(input->readX(0).size(), output->readX(0).size());
+
+ TS_ASSERT(output->run().hasProperty("NumAllSpectra"))
+ TS_ASSERT(output->run().hasProperty("NumMaskSpectra"))
+ TS_ASSERT(output->run().hasProperty("NumZeroSpectra"))
+
}
void testRebinnedOutputSum()
@@ -195,6 +222,7 @@ public:
std::string outName = "rebin_sum";
AnalysisDataService::Instance().addOrReplace(inName, ws);
+ size_t nHist = ws->getNumberHistograms();
// Start with a clean algorithm
Mantid::Algorithms::SumSpectra alg3;
@@ -223,9 +251,91 @@ public:
TS_ASSERT_DELTA(output->dataY(0)[5], 0.66666, 1.e-5);
TS_ASSERT_DELTA(output->dataE(0)[5], 0.47140452079103173, 1.e-5);
TS_ASSERT_EQUALS(output->dataF(0)[5], 3.);
+
+ TS_ASSERT(output->run().hasProperty("NumAllSpectra"))
+ TS_ASSERT(output->run().hasProperty("NumMaskSpectra"))
+ TS_ASSERT(output->run().hasProperty("NumZeroSpectra"))
+
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(nHist),output->run().getLogData("NumAllSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(0),output->run().getLogData("NumMaskSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(3),output->run().getLogData("NumZeroSpectra")->value())
+
}
+ void testExecNoLimitsWeighted()
+ {
+ Mantid::Algorithms::SumSpectra alg2;
+ TS_ASSERT_THROWS_NOTHING( alg2.initialize());
+ TS_ASSERT( alg2.isInitialized() );
+
+
+ const Mantid::MantidVec &y0 = inputSpace->readY(0);
+ const Mantid::MantidVec &e0 = inputSpace->readE(0);
+ // Set the properties
+ alg2.setProperty("InputWorkspace",inputSpace);
+ const std::string outputSpace2 = "SumSpectraOut2";
+ alg2.setPropertyValue("OutputWorkspace",outputSpace2);
+ alg2.setProperty("IncludeMonitors",false);
+ alg2.setProperty("WeightedSum",true);
+
+ TS_ASSERT_THROWS_NOTHING( alg2.execute());
+ TS_ASSERT( alg2.isExecuted() );
+
+ // Get back the saved workspace
+ Workspace_sptr output;
+ TS_ASSERT_THROWS_NOTHING(output = AnalysisDataService::Instance().retrieve(outputSpace2));
+ Workspace2D_const_sptr output2D = boost::dynamic_pointer_cast<const Workspace2D>(output);
+
+ TS_ASSERT_EQUALS( output2D->getNumberHistograms(), 1);
+
+ const Mantid::MantidVec &x = output2D->readX(0);
+ const Mantid::MantidVec &y = output2D->readY(0);
+ const Mantid::MantidVec &e = output2D->readE(0);
+ TS_ASSERT_EQUALS( x.size(), 103 );
+ TS_ASSERT_EQUALS( y.size(), 102 );
+ TS_ASSERT_EQUALS( e.size(), 102 );
+
+ TS_ASSERT(output2D->run().hasProperty("NumAllSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumMaskSpectra"))
+ TS_ASSERT(output2D->run().hasProperty("NumZeroSpectra"))
+
+ // Spectra at workspace index 1 is masked, 8 & 9 are monitors
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(nTestHist-3),output2D->run().getLogData("NumAllSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(1),output2D->run().getLogData("NumMaskSpectra")->value())
+ TS_ASSERT_EQUALS(boost::lexical_cast<std::string>(0),output2D->run().getLogData("NumZeroSpectra")->value())
+
+ size_t nSignals = nTestHist-3;
+ // Check a few bins
+ TS_ASSERT_EQUALS( x[0], inputSpace->readX(0)[0] );
+ TS_ASSERT_EQUALS( x[50], inputSpace->readX(0)[50] );
+ TS_ASSERT_EQUALS( x[100], inputSpace->readX(0)[100] );
+ TS_ASSERT_EQUALS( y[7], y0[7] );
+ TS_ASSERT_EQUALS( y[38],y0[38] );
+ TS_ASSERT_EQUALS( y[72],y0[72] );
+ TS_ASSERT_DELTA( e[28], std::sqrt(double(nSignals))*e0[28], 0.00001 );
+ TS_ASSERT_DELTA( e[47], std::sqrt(double(nSignals))*e0[47], 0.00001 );
+ TS_ASSERT_DELTA( e[99], std::sqrt(double(nSignals))*e0[99], 0.00001 );
+
+
+ // Check the detectors mapped to the single spectra
+ const ISpectrum * spec = output2D->getSpectrum(0);
+ TS_ASSERT_EQUALS( spec->getSpectrumNo(), 1);
+ // Spectra at workspace index 1 is masked, 8 & 9 are monitors
+ TS_ASSERT_EQUALS( spec->getDetectorIDs().size(), 7);
+ TS_ASSERT( spec->hasDetectorID(1) );
+ TS_ASSERT( spec->hasDetectorID(3) );
+ TS_ASSERT( spec->hasDetectorID(4) );
+ TS_ASSERT( spec->hasDetectorID(5) );
+ TS_ASSERT( spec->hasDetectorID(6) );
+ TS_ASSERT( spec->hasDetectorID(7) );
+
+
+
+
+ }
+
private:
+ int nTestHist;
Mantid::Algorithms::SumSpectra alg; // Test with range limits
MatrixWorkspace_sptr inputSpace;
};
--
GitLab