Skip to content

Commit

Permalink
Refs #4720: NormaliseToMonitor preserves events
Browse files Browse the repository at this point in the history
  • Loading branch information
Janik Zikovsky committed Apr 4, 2012
1 parent cb40aa8 commit 6fc6962
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 31 deletions.
85 changes: 60 additions & 25 deletions Code/Mantid/Framework/Algorithms/src/NormaliseToMonitor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,9 @@ In both cases, the [[Divide]] algorithm is used to perform the normalisation.

#include <cfloat>
#include <iomanip>
#include "MantidDataObjects/EventWorkspace.h"

using namespace Mantid::DataObjects;

namespace Mantid
{
Expand Down Expand Up @@ -426,11 +429,16 @@ API::MatrixWorkspace_sptr NormaliseToMonitor::extractMonitorSpectrum(API::Matrix
IAlgorithm_sptr childAlg = createSubAlgorithm("ExtractSingleSpectrum");
childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
childAlg->setProperty<int>("WorkspaceIndex", static_cast<int>(index));

childAlg->executeAsSubAlg();
MatrixWorkspace_sptr outWS = childAlg->getProperty("OutputWorkspace");

IAlgorithm_sptr alg = createSubAlgorithm("ConvertToMatrixWorkspace");
alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outWS);
alg->executeAsSubAlg();
outWS = alg->getProperty("OutputWorkspace");

// Only get to here if successful
return childAlg->getProperty("OutputWorkspace");
return outWS;
}

/** Sets the maximum and minimum X values of the monitor spectrum to use for integration
Expand Down Expand Up @@ -513,9 +521,25 @@ void NormaliseToMonitor::normaliseByIntegratedCount(API::MatrixWorkspace_sptr in
void NormaliseToMonitor::normaliseBinByBin(API::MatrixWorkspace_sptr inputWorkspace,
API::MatrixWorkspace_sptr& outputWorkspace)
{
EventWorkspace_sptr inputEvent = boost::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
EventWorkspace_sptr outputEvent;

// Only create output workspace if different to input one
if (outputWorkspace != inputWorkspace )
outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
if (outputWorkspace != inputWorkspace )
{
if (inputEvent)
{
//Make a brand new EventWorkspace
outputEvent = boost::dynamic_pointer_cast<EventWorkspace>(
API::WorkspaceFactory::Instance().create("EventWorkspace", inputEvent->getNumberHistograms(), 2, 1));
//Copy geometry and data
API::WorkspaceFactory::Instance().initializeFromParent(inputEvent, outputEvent, false);
outputEvent->copyDataFrom( (*inputEvent) );
outputWorkspace = boost::dynamic_pointer_cast<MatrixWorkspace>(outputEvent);
}
else
outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
}

// Get hold of the monitor spectrum
const MantidVec& monX = m_monitor->readX(0);
Expand All @@ -524,6 +548,7 @@ void NormaliseToMonitor::normaliseBinByBin(API::MatrixWorkspace_sptr inputWorksp
// Calculate the overall normalisation just the once if bins are all matching
if (m_commonBins) this->normalisationFactor(m_monitor->readX(0),&monY,&monE);


const size_t numHists = inputWorkspace->getNumberHistograms();
MantidVec::size_type specLength = inputWorkspace->blocksize();
Progress prog(this,0.0,1.0,numHists);
Expand All @@ -549,31 +574,41 @@ void NormaliseToMonitor::normaliseBinByBin(API::MatrixWorkspace_sptr inputWorksp
this->normalisationFactor(X,Y,E);
}

const MantidVec& inY = inputWorkspace->readY(i);
const MantidVec& inE = inputWorkspace->readE(i);
MantidVec& YOut = outputWorkspace->dataY(i);
MantidVec& EOut = outputWorkspace->dataE(i);
outputWorkspace->dataX(i) = inputWorkspace->readX(i);
// The code below comes more or less straight out of Divide.cpp
for (MantidVec::size_type k = 0; k < specLength; ++k)
if (inputEvent)
{
// ----------------------------------- EventWorkspace ---------------------------------------
EventList & outEL = outputEvent->getEventList(i);
outEL.divide(X, *Y, *E);
}
else
{
// Get references to the input Y's
const double& leftY = inY[k];
const double& rightY = (*Y)[k];
// ----------------------------------- Workspace2D ---------------------------------------
const MantidVec& inY = inputWorkspace->readY(i);
const MantidVec& inE = inputWorkspace->readE(i);
MantidVec& YOut = outputWorkspace->dataY(i);
MantidVec& EOut = outputWorkspace->dataE(i);
outputWorkspace->dataX(i) = inputWorkspace->readX(i);
// The code below comes more or less straight out of Divide.cpp
for (MantidVec::size_type k = 0; k < specLength; ++k)
{
// Get references to the input Y's
const double& leftY = inY[k];
const double& rightY = (*Y)[k];

// Calculate result and store in local variable to avoid overwriting original data if
// output workspace is same as one of the input ones
const double newY = leftY/rightY;
// Calculate result and store in local variable to avoid overwriting original data if
// output workspace is same as one of the input ones
const double newY = leftY/rightY;

if (fabs(rightY)>1.0e-12 && fabs(newY)>1.0e-12)
{
const double lhsFactor = (inE[k]<1.0e-12|| fabs(leftY)<1.0e-12) ? 0.0 : pow((inE[k]/leftY),2);
const double rhsFactor = (*E)[k]<1.0e-12 ? 0.0 : pow(((*E)[k]/rightY),2);
EOut[k] = std::abs(newY) * sqrt(lhsFactor+rhsFactor);
}
if (fabs(rightY)>1.0e-12 && fabs(newY)>1.0e-12)
{
const double lhsFactor = (inE[k]<1.0e-12|| fabs(leftY)<1.0e-12) ? 0.0 : pow((inE[k]/leftY),2);
const double rhsFactor = (*E)[k]<1.0e-12 ? 0.0 : pow(((*E)[k]/rightY),2);
EOut[k] = std::abs(newY) * sqrt(lhsFactor+rhsFactor);
}

// Now store the result
YOut[k] = newY;
// Now store the result
YOut[k] = newY;
} // end Workspace2D case
} // end loop over current spectrum

if (!m_commonBins) { delete Y; delete E; }
Expand Down
37 changes: 31 additions & 6 deletions Code/Mantid/Framework/Algorithms/test/NormaliseToMonitorTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#include "MantidAPI/SpectraDetectorMap.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/Property.h"
#include "MantidKernel/SingletonHolder.h"
#include "MantidAPI/FrameworkManager.h"

using namespace Mantid::Kernel;
using namespace Mantid::API;
Expand All @@ -30,14 +32,18 @@ class NormaliseToMonitorTestHelper : public NormaliseToMonitor
class NormaliseToMonitorTest : public CxxTest::TestSuite
{
private:
NormaliseToMonitorTestHelper norm;


public:

static NormaliseToMonitorTest *createSuite() { return new NormaliseToMonitorTest(); }
static void destroySuite(NormaliseToMonitorTest *suite) { delete suite; }

NormaliseToMonitorTest()
{
}

void setUp()
{
MatrixWorkspace_sptr input = WorkspaceCreationHelper::Create2DWorkspace123(3,10,1);
// Change the data in the monitor spectrum
Expand Down Expand Up @@ -66,7 +72,7 @@ class NormaliseToMonitorTest : public CxxTest::TestSuite
instr->markAsDetector(det);
input->replaceSpectraMap(new SpectraDetectorMap(forSpecDetMap, forSpecDetMap, 3));

AnalysisDataService::Instance().add("normMon",input);
AnalysisDataService::Instance().addOrReplace("normMon",input);

// Create a single spectrum workspace to be the monitor one
MatrixWorkspace_sptr monWS = WorkspaceCreationHelper::Create2DWorkspaceBinned(1,20,0.1,0.5);
Expand All @@ -82,33 +88,42 @@ class NormaliseToMonitorTest : public CxxTest::TestSuite
//instr->markAsMonitor(mon2);
monWS->replaceSpectraMap(new SpectraDetectorMap(forSpecDetMap2, forSpecDetMap2, 1));

AnalysisDataService::Instance().add("monWS",monWS);
AnalysisDataService::Instance().addOrReplace("monWS",monWS);
}

void testName()
{
NormaliseToMonitorTestHelper norm;
TS_ASSERT_EQUALS( norm.name(), "NormaliseToMonitor" )
}

void testVersion()
{
NormaliseToMonitorTestHelper norm;
TS_ASSERT_EQUALS( norm.version(), 1 )
}

void testInit()
{
NormaliseToMonitorTestHelper norm;
TS_ASSERT_THROWS_NOTHING( norm.initialize() )
TS_ASSERT( norm.isInitialized() )
}

void testExec()
void dotestExec(bool events)
{
if ( !norm.isInitialized() ) norm.initialize();
NormaliseToMonitorTestHelper norm;
if (events)
FrameworkManager::Instance().exec("ConvertToEventWorkspace", 8,
"InputWorkspace", "normMon",
"GenerateZeros", "1",
"GenerateMultipleEvents", "0",
"OutputWorkspace", "normMon");

if ( !norm.isInitialized() ) norm.initialize();
// Check it fails if properties haven't been set
TS_ASSERT_THROWS( norm.execute(), std::runtime_error )
TS_ASSERT( ! norm.isExecuted() )

TS_ASSERT_THROWS_NOTHING( norm.setPropertyValue("InputWorkspace","normMon") )
TS_ASSERT_THROWS_NOTHING( norm.setPropertyValue("OutputWorkspace","normMon2") )
TS_ASSERT_THROWS_NOTHING( norm.setPropertyValue("MonitorSpectrum","0") )
Expand Down Expand Up @@ -139,6 +154,16 @@ class NormaliseToMonitorTest : public CxxTest::TestSuite
}


void testExec()
{
dotestExec(false);
}

void testExec_Events()
{
dotestExec(true);
}


void testNormaliseByIntegratedCount()
{
Expand Down

0 comments on commit 6fc6962

Please sign in to comment.