Skip to content

Commit

Permalink
Refs #4814. FindeDetectorsOutsideLimits outputs a SpecialWorkspace2D.
Browse files Browse the repository at this point in the history
  • Loading branch information
peterfpeterson committed Feb 21, 2012
1 parent 66cef25 commit ebd70bc
Show file tree
Hide file tree
Showing 4 changed files with 53 additions and 35 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,8 @@ namespace Mantid
const int indexMax, const double lower,
const double upper, const bool outputWorkspace2D = false);

API::MatrixWorkspace_sptr generateEmptyMask(API::MatrixWorkspace_sptr inputWS, const bool initialize = true);

/// Calculate the median of the given workspace. This assumes that the input workspace contains
/// integrated counts
double calculateMedian(API::MatrixWorkspace_sptr input, bool excludeZeroes);
Expand Down
37 changes: 36 additions & 1 deletion Code/Mantid/Framework/Algorithms/src/DetectorDiagnostic.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Exception.h"
#include "MantidDataObjects/EventWorkspaceHelpers.h"
#include "MantidDataObjects/SpecialWorkspace2D.h"
#include <boost/math/special_functions/fpclassify.hpp>
#include <gsl/gsl_statistics.h>
#include <cfloat>
Expand Down Expand Up @@ -80,7 +81,41 @@ namespace Mantid
return finalOutputW;
}


/**
* Create a masking workspace to return.
*
* @param inputWS The workspace to initialize from. The instrument is copied from this.
* @param initialize Whether or not to set all of the values to keep the data.
*/
MatrixWorkspace_sptr DetectorDiagnostic::generateEmptyMask(API::MatrixWorkspace_sptr inputWS, const bool initialize)
{
// Create a new workspace for the results, copy from the input to ensure that we copy over the instrument and current masking
DataObjects::SpecialWorkspace2D* maskWS = new DataObjects::SpecialWorkspace2D();
maskWS->initialize(inputWS->getNumberHistograms(), 1, 1);
MatrixWorkspace_sptr outputWS(maskWS);
WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, false);
outputWS->setTitle(inputWS->getTitle());

if (initialize) // initialize the mask as keep everything
{
// initialize the mask as keep everything
const double liveValue(0.0); // keep the data - delete is 1.0
const double errors(0.0); // uncertainty
int64_t size = outputWS->getNumberHistograms();
PARALLEL_FOR1(outputWS)
for (int64_t i = 0; i < size; i++)
{
PARALLEL_START_INTERUPT_REGION
outputWS->dataY(i)[0] = liveValue; // by default it stays
outputWS->dataE(i)[0] = errors;
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}

return outputWS;
}

/**
* Fnds the median of values in single bin histograms rejecting spectra from masked
* detectors and the results of divide by zero (infinite and NaN).
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -118,18 +118,16 @@ namespace Mantid
if (boost::dynamic_pointer_cast<EventWorkspace>(countsWS))
throw std::runtime_error("Error! Integration output is not a Workspace2D.");

// This will become the output workspace
countsWS->isDistribution(false);
countsWS->setYUnit("");
countsWS->instrumentParameters();
// Create a new workspace for the results, copy from the input to ensure that we copy over the instrument and current masking
MatrixWorkspace_sptr outputWS = this->generateEmptyMask(inputWS);

const double liveValue(1.0);
const double deadValue(1.0); // delete the data

const int diagLength = static_cast<int>(countsWS->getNumberHistograms());
const int progStep = static_cast<int>(std::ceil(diagLength / 100.0));

int numFailed(0);
PARALLEL_FOR1(countsWS)
PARALLEL_FOR2(countsWS, outputWS)
for (int i = 0; i < diagLength; ++i)
{
if (i % progStep == 0)
Expand All @@ -148,7 +146,7 @@ namespace Mantid
// Mark no detector spectra as failed
if( !det )
{
countsWS->maskWorkspaceIndex(i);
outputWS->dataY(i)[0] = deadValue;
PARALLEL_ATOMIC
++numFailed;
continue;
Expand All @@ -157,37 +155,32 @@ namespace Mantid
if( det->isMonitor() )
{
// Don't include but don't mask either
countsWS->dataY(i)[0] = liveValue;
continue;
}

const double & yValue = countsWS->readY(i)[0];
// Mask out NaN and infinite
if( boost::math::isinf(yValue) || boost::math::isnan(yValue) )
{
countsWS->maskWorkspaceIndex(i);
outputWS->dataY(i)[0] = deadValue;
PARALLEL_ATOMIC
++numFailed;
continue;
}

if ( yValue <= lowThreshold || yValue >= highThreshold)
{
countsWS->maskWorkspaceIndex(i);
outputWS->dataY(i)[0] = deadValue;
PARALLEL_ATOMIC
++numFailed;
}
else
{
// Value passed the tests, flag it.
countsWS->dataY(i)[0] = liveValue;
}

}

g_log.information() << numFailed << " spectra fell outside the given limits.\n";
setProperty("NumberOfFailures", numFailed);
// Assign it to the output workspace property
setProperty("OutputWorkspace", countsWS);
setProperty("OutputWorkspace", outputWS);
}

}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -93,25 +93,19 @@ class FindDetectorsOutsideLimitsTest : public CxxTest::TestSuite
const int numFailed = alg.getProperty("NumberOfFailures");
TS_ASSERT_EQUALS(numFailed, 11);

const double liveValue(1.0);
const double maskValue(0.0);
const double liveValue(0.0);
const double maskValue(1.0);
for (int i=0; i< sizey; i++)
{
const double val = work_out->readY(i)[0];
double valExpected = liveValue;
// Check masking
IDetector_const_sptr det;
TS_ASSERT_THROWS_NOTHING(det = work_out->getDetector(i));
bool maskExpected(false);
// Spectra set up with yVeryDead fail low counts or yStrange fail on high
if ( i%2 == 0 || i == 19 )
{
valExpected = maskValue;
maskExpected = true;
}
if(det)
{
TS_ASSERT_EQUALS(det->isMasked(), maskExpected);
}

TS_ASSERT_DELTA(val,valExpected,1e-9);
Expand All @@ -136,16 +130,10 @@ class FindDetectorsOutsideLimitsTest : public CxxTest::TestSuite
// Check masking
IDetector_const_sptr det;
TS_ASSERT_THROWS_NOTHING(det = work_out->getDetector(i));
bool maskExpected(false);
// Spectra set up with yVeryDead fail low counts or yStrange fail on high
if ( i%2 == 0 )
{
valExpected = maskValue;
maskExpected = true;
}
if(det)
{
TS_ASSERT_EQUALS(det->isMasked(), maskExpected);
}

TS_ASSERT_DELTA(val,valExpected,1e-9);
Expand Down Expand Up @@ -183,10 +171,10 @@ class FindDetectorsOutsideLimitsTest : public CxxTest::TestSuite
MatrixWorkspace_sptr work_out;
TS_ASSERT_THROWS_NOTHING(work_out = boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve("testdead_out")));

TS_ASSERT_EQUALS( work_out->dataY(0)[0], 1.0);
TS_ASSERT_EQUALS( work_out->dataY(9)[0], 1.0);
TS_ASSERT_EQUALS( work_out->dataY(10)[0], 0.0);
TS_ASSERT_EQUALS( work_out->dataY(11)[0], 1.0);
TS_ASSERT_EQUALS( work_out->dataY(0)[0], 0.0);
TS_ASSERT_EQUALS( work_out->dataY(9)[0], 0.0);
TS_ASSERT_EQUALS( work_out->dataY(10)[0], 1.0);
TS_ASSERT_EQUALS( work_out->dataY(11)[0], 0.0);

AnalysisDataService::Instance().remove("testdead_in");
AnalysisDataService::Instance().remove("testdead_out");
Expand Down

0 comments on commit ebd70bc

Please sign in to comment.