Skip to content

Commit

Permalink
Implementing generating masking workspace. Refs #4370.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Dec 29, 2011
1 parent 9f5f764 commit c2e67af
Show file tree
Hide file tree
Showing 5 changed files with 269 additions and 87 deletions.
1 change: 1 addition & 0 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FindDetectorsOutsideLimits.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/DetectorDiagnostic.h"
#include "MantidDataObjects/SpecialWorkspace2D.h"

namespace Mantid
{
Expand Down
7 changes: 4 additions & 3 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MedianDetectorTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/DetectorDiagnostic.h"
#include "MantidDataObjects/SpecialWorkspace2D.h"

namespace Mantid
{
Expand Down Expand Up @@ -83,9 +84,9 @@ namespace Mantid
API::MatrixWorkspace_sptr getSolidAngles(int firstSpec, int lastSpec);
/// Do the tests and mask those that fail
int doDetectorTests(const API::MatrixWorkspace_sptr countWorkspace,
API::MatrixWorkspace_sptr maskWS,
const double average,
const std::set<int> & badIndices);
DataObjects::SpecialWorkspace2D_sptr maskWS,
const double average,
const std::set<int> & badIndices);

/// Input workspace
API::MatrixWorkspace_sptr m_inputWS;
Expand Down
78 changes: 59 additions & 19 deletions Code/Mantid/Framework/Algorithms/src/FindDetectorsOutsideLimits.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,14 @@ namespace Mantid
new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input),
"Name of the input workspace2D" );
declareProperty(
new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output),
new WorkspaceProperty<DataObjects::SpecialWorkspace2D>("OutputWorkspace","",Direction::Output),
"Each histogram from the input workspace maps to a histogram in this\n"
"workspace with one value that indicates if there was a dead detector" );
"SpecialWorkspace2D (Masking workdpace) with one value that indicates if there was a dead detector" );
/*
declareProperty(
new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace2", "", Direction::Output),
"Same as OutputWorksapce, but in Workspace2D");
*/
declareProperty("HighThreshold", 0.0,
"Spectra whose total number of counts are equal to or above this value\n"
"will be marked bad (default 0)" );
Expand All @@ -71,6 +76,7 @@ namespace Mantid
*/
void FindDetectorsOutsideLimits::exec()
{
// 1. Get properties (inputs)
double lowThreshold = getProperty("LowThreshold");
double highThreshold = getProperty("HighThreshold");
if ( highThreshold <= lowThreshold )
Expand All @@ -85,53 +91,84 @@ namespace Mantid
const double rangeLower = getProperty("RangeLower");
const double rangeUpper = getProperty("RangeUpper");

// Get the integrated input workspace; converting to a Workspace2D
// 2. Create output workspace and set it up & check input workspace size
Geometry::Instrument_const_sptr instrument = inputWS->getInstrument();
DataObjects::SpecialWorkspace2D_sptr maskWS(new DataObjects::SpecialWorkspace2D(instrument));
this->setProperty("OutputWorkspace", maskWS);

if (maskWS->getNumberHistograms() > inputWS->getNumberHistograms()){
g_log.warning() << "Input Workspace " << inputWS->getName() << " has different number of spectra (" <<
inputWS->getNumberHistograms() << ") to new generated masking workspace (" <<
maskWS->getNumberHistograms() << "). " << " Input workspace might be focused" << std::endl;
}
else if (maskWS->getNumberHistograms() < inputWS->getNumberHistograms()){
g_log.error() << "Input Workspace " << inputWS->getName() << " has different number of spectra (" <<
inputWS->getNumberHistograms() << ") to new generated masking workspace (" <<
maskWS->getNumberHistograms() << "). " << " Input workspace might be focused" << std::endl;
throw std::runtime_error("Error! This cannot be right. Instrument must be wrong.");
}

// 3. Get the integrated input workspace; converting to a Workspace2D
MatrixWorkspace_sptr integratedWorkspace =
integrateSpectra(inputWS, 0, EMPTY_INT(), rangeLower, rangeUpper, true);

EventWorkspace_sptr ew = boost::dynamic_pointer_cast<EventWorkspace>( integratedWorkspace);
EventWorkspace_sptr ew = boost::dynamic_pointer_cast<EventWorkspace>(integratedWorkspace);
if (ew)
throw std::runtime_error("Error! Integration output is not a Workspace2D.");

//set the workspace to have no units
integratedWorkspace->isDistribution(false);
integratedWorkspace->setYUnit("");

// i. set the workspace to have no units
// TODO remove the setup on integratedWorkspace after testing
// integratedWorkspace->isDistribution(false);
// integratedWorkspace->setYUnit("");

// 4. Read integrated and set up the Masking workspace
const double liveValue(1.0);
// iterate over the data values setting the live and dead values
const double deadValue = 0.0;

// a) iterate over the data values setting the live and dead values
const int numSpec = static_cast<int>(integratedWorkspace->getNumberHistograms());
const int progStep = static_cast<int>(ceil(numSpec / 100.0));

// Keep track of those reading low/high
// b) Keep track of those reading low/high
int numLow(0), numHigh(0);

// c) Loop
for (int i = 0; i < numSpec; ++i)
{
// i) Progress
if (i % progStep == 0)
{
progress(static_cast<double>(i)/numSpec);
interruption_point();
}

// ii) Set up value
const double & yValue = integratedWorkspace->readY(i)[0];
if ( yValue <= lowThreshold )
{
// Value failed
integratedWorkspace->maskWorkspaceIndex(i);

maskWS->dataY(i)[0] = deadValue;

++numLow;
continue;
}

if ( yValue >= highThreshold )
} else if ( yValue >= highThreshold )
{
// Value failed
integratedWorkspace->maskWorkspaceIndex(i);

maskWS->dataY(i)[0] = deadValue;

++numHigh;
continue;
}

// Value passed the tests, flag it.
integratedWorkspace->dataY(i)[0] = liveValue;
}
} else {
// Value passed the tests, flag it.
integratedWorkspace->dataY(i)[0] = liveValue;

maskWS->dataY(i)[0] = liveValue;
} // ENDIFELSE
} // ENDFOR

const int totalNumFailed(numLow+numHigh);
g_log.notice() << "Found a total of " << totalNumFailed << " failing "
Expand All @@ -140,7 +177,10 @@ namespace Mantid

setProperty("NumberOfFailures", totalNumFailed);
// Assign it to the output workspace property
setProperty("OutputWorkspace", integratedWorkspace);
// setProperty("OutputWorkspace2", integratedWorkspace);
setProperty("OutputWorkspace", maskWS);

return;
}

}
Expand Down
130 changes: 65 additions & 65 deletions Code/Mantid/Framework/Algorithms/src/MedianDetectorTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ namespace Mantid
new HistogramValidator<>),
"Name of the input workspace" );
declareProperty(
new WorkspaceProperty<>("OutputWorkspace","",Direction::Output),
new WorkspaceProperty<DataObjects::SpecialWorkspace2D>("OutputWorkspace","",Direction::Output),
"A MaskWorkspace where 0 denotes a masked spectra. Any spectra containing"
"a zero is also masked on the output");

Expand Down Expand Up @@ -90,43 +90,39 @@ namespace Mantid
*/
void MedianDetectorTest::exec()
{
// 1. Get input arguments
retrieveProperties();



//Adds the counts from all the bins and puts them in one total bin
// 2. Adds the counts from all the bins and puts them in one total bin
MatrixWorkspace_sptr counts = integrateSpectra(m_inputWS, m_minSpec, m_maxSpec,
m_rangeLower, m_rangeUpper, true);
counts = convertToRate(counts);

// FIXME: The next section that calculates the solid angle is commented out until
// the SolidAngle algorithm is corrected to return the correct number.
// (see http://trac.mantidproject.org/mantid/ticket/2596)
// FIXME: The next section that calculates the solid angle is commented out until
// the SolidAngle algorithm is corrected to return the correct number.
// (see http://trac.mantidproject.org/mantid/ticket/2596)

//MatrixWorkspace_sptr angles = getSolidAngles(m_minSpec, m_maxSpec);

//Gets the count rate per solid angle (in steradians), if it exists, for each spectrum
//this calculation is optional, it depends on angle information existing
//if ( angles.use_count() == 1 )
//{
//if some numbers in angles are zero we will get the infinity flag value
// in the output work space which needs to be dealt with later
//counts = counts/angles;
//if some numbers in angles are zero we will get the infinity flag value
// in the output work space which needs to be dealt with later
//counts = counts/angles;
//}

// End of Solid Angle commented out section


// An average of the data, the median is less influenced by a small number of huge values than the mean
// 3. An average of the data, the median is less influenced by a small number of huge values than the mean
std::set<int> badIndices;
double average = calculateMedian(counts, badIndices);

// Create a workspace for the output
MatrixWorkspace_sptr maskWS = WorkspaceFactory::Instance().create(counts);
// Make sure the output is simple
maskWS->isDistribution(false);
maskWS->setYUnit("");
maskWS->instrumentParameters().clear();
// 4. Initialize output masking workspace
Geometry::Instrument_const_sptr myinstrument = m_inputWS->getInstrument();
DataObjects::SpecialWorkspace2D_sptr maskWS(new DataObjects::SpecialWorkspace2D(myinstrument));

int numFailed = doDetectorTests(counts, maskWS, average, badIndices);

Expand Down Expand Up @@ -227,8 +223,8 @@ namespace Mantid
* @return The number of detectors that failed the tests, not including those skipped
*/
int MedianDetectorTest::doDetectorTests(const API::MatrixWorkspace_sptr countWorkspace,
API::MatrixWorkspace_sptr maskWS,
const double average, const std::set<int> & badIndices)
DataObjects::SpecialWorkspace2D_sptr maskWS,
const double average, const std::set<int> & badIndices)
{
g_log.information("Applying the criteria to find failing detectors");

Expand All @@ -242,83 +238,87 @@ namespace Mantid
const int numSpec(m_maxSpec - m_minSpec);
const int progStep = static_cast<int>(ceil(numSpec/30.0));

const double live_value(1.0);
int numLow(0), numHigh(0);
size_t numBadPixels = 0;

//set the workspace to have no units
countWorkspace->isDistribution(false);
countWorkspace->setYUnit("");

double liveValue = 1.0;
double deadValue = 0.0;

PARALLEL_FOR1(countWorkspace)
for (int i = 0; i <= numSpec; ++i)
{
PARALLEL_START_INTERUPT_REGION

// update the progressbar information
// i. update the progressbar information
if (i % progStep == 0)
{
progress(advanceProgress(progStep*static_cast<double>(RTMarkDetects)/numSpec));
}


// If the value is not in the badIndices set then assume it is good
// else skip tests for it
// ii. If the value is not in the badIndices set then assume it is good
// else skip tests for it
if( badIndices.count(i) == 1 )
{
maskWS->maskWorkspaceIndex(i);
continue;
}
//Do tests
const double sig = minSigma*countWorkspace->readE(i)[0];
// Check the significance value is okay
if( boost::math::isinf(std::abs(sig)) || boost::math::isinf(sig) )
{
PARALLEL_CRITICAL(MedianDetectorTest_failed_a)
{
maskWS->maskWorkspaceIndex(i);
++numLow;
}
continue;
}
// Bad pixel
maskWS->dataY(i)[0] = deadValue;
numBadPixels ++;

const double yIn = countWorkspace->dataY(i)[0];
if ( yIn <= lowLim )
{
// compare the difference against the size of the errorbar -statistical significance check
if(average - yIn > sig)
} else {
//Good pixel: do tests

const double sig = minSigma*countWorkspace->readE(i)[0];

// (a) Check the significance value is okay
if( boost::math::isinf(std::abs(sig)) || boost::math::isinf(sig) )
{
PARALLEL_CRITICAL(MedianDetectorTest_failed_b)
// sig is infinity
PARALLEL_CRITICAL(MedianDetectorTest_failed_a)
{
maskWS->maskWorkspaceIndex(i);
maskWS->dataY(i)[0] = deadValue;
++numLow;
}
continue;
}
}
if (yIn >= highLim)
{
// compare the difference against the size of the errorbar -statistical significance check
if(yIn - average > sig)
{
PARALLEL_CRITICAL(MedianDetectorTest_failed_c)
} else {
// sig is not infinity
const double yIn = countWorkspace->dataY(i)[0];
if ( (yIn <= lowLim) && (average - yIn > sig))
{
maskWS->maskWorkspaceIndex(i);
++numHigh;
// compare the difference against the size of the errorbar -statistical significance check
PARALLEL_CRITICAL(MedianDetectorTest_failed_b)
{
maskWS->dataY(i)[0] = deadValue;
++numLow;
}
}
continue;
}
}
// Reaching here passes the tests
maskWS->dataY(i)[0] = live_value;
else if ( (yIn >= highLim) && (yIn - average > sig) )
{
// compare the difference against the size of the errorbar -statistical significance check
PARALLEL_CRITICAL(MedianDetectorTest_failed_c)
{
maskWS->dataY(i)[0] = deadValue;
++numHigh;
}
} else {
// Reaching here passes the tests
maskWS->dataY(i)[0] = liveValue;

} // IF-ELSE: test
} // IF-ELSE: Sig is infinity?
} // IF-ELSE: Bad/good pixels

PARALLEL_END_INTERUPT_REGION
}
} // ENDFOR
PARALLEL_CHECK_INTERUPT_REGION

// Log finds
g_log.information() << "-- Detector tests --\n "
<< "Number recording low: " << numLow << "\n"
<< "Number recording high: " << numHigh << "\n";
<< "Number recording high: " << numHigh << "\n"
<< "Number bad pixels: << " << numBadPixels << std::endl;

return (numLow + numHigh);
}
Expand Down

0 comments on commit c2e67af

Please sign in to comment.