Skip to content

Commit

Permalink
Refs #9095. Enabled to filter out peaks by resolution.
Browse files Browse the repository at this point in the history 
Also added unit tests of this new feature (not passed yet)
  • Loading branch information
@wdzhou
wdzhou committed May 20, 2014
1 parent ca8bd18 commit b968c07
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 4 deletions.
5 changes: 4 additions & 1 deletion Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,10 @@ class DLLExport GetDetOffsetsMultiPeaks: public API::Algorithm
API::MatrixWorkspace_const_sptr m_inputResolutionWS;
/// Flag of use input resolution WS
bool m_hasInputResolution;
/// Lower boundary of allowed peak width as resolution
double m_minResFactor;
/// Upper boundary of allowed peak width as resolution
double m_maxResFactor;

DataObjects::OffsetsWorkspace_sptr outputW;
/// Output workspace for debugging purpose
Expand All @@ -178,7 +182,6 @@ class DLLExport GetDetOffsetsMultiPeaks: public API::Algorithm
std::vector<std::vector<double> > m_vecFitWindow;



};

} // namespace Algorithm
Expand Down
34 changes: 31 additions & 3 deletions Code/Mantid/Framework/Algorithms/src/GetDetOffsetsMultiPeaks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,10 @@ namespace Algorithms
auto ddodwsprop = new WorkspaceProperty<MatrixWorkspace>("FittedResolutionWorkspace", "ResolutionWS", Direction::Output);
declareProperty(ddodwsprop, "Name of the resolution workspace containing delta(d)/d for each unmasked spectrum. ");

declareProperty("MinimumResolutionFactor", 0.1, "Factor of the minimum allowed Delta(d)/d of any peak to its suggested Delta(d)/d. ");

declareProperty("MiximumResolutionFactor", 10.0, "Factor of the maximum allowed Delta(d)/d of any peak to its suggested Delta(d)/d. ");

}


Expand Down Expand Up @@ -443,6 +447,12 @@ namespace Algorithms
m_inputResolutionWS = getProperty("InputResolutionWorkspace");
m_hasInputResolution = true;

m_minResFactor = getProperty("MinimumResolutionFactor");
m_maxResFactor = getProperty("MaximumResolutionFactor");

if (m_minResFactor >= m_maxResFactor)
throw std::runtime_error("Input peak resolution boundary is 0 or negative.");

// Check
if (m_inputResolutionWS->getNumberHistograms() != m_inputWS->getNumberHistograms())
throw std::runtime_error("Input workspace does not match resolution workspace. ");
Expand Down Expand Up @@ -564,11 +574,13 @@ namespace Algorithms
maskWS->dataY(workspaceIndex)[0] = offsetresult.mask;

// check the average value of delta(d)/d. if it is far off the theorical value, output
// FIXME - This warning should not appear by filtering out peaks that are too wide or narrow.
// TODO - Delete the if statement below if it is never triggered.
if (m_hasInputResolution)
{
double pixelresolution = m_inputResolutionWS->readY(wi)[0];
if (offsetresult.resolution > 10*pixelresolution || offsetresult.resolution < 0.1*pixelresolution)
g_log.notice() << "Spectrum " << wi << " delta(d)/d = " << offsetresult.resolution << "\n";
g_log.warning() << "Spectrum " << wi << " delta(d)/d = " << offsetresult.resolution << "\n";
}
}
} // ENDFOR (detectors)
Expand Down Expand Up @@ -1074,6 +1086,23 @@ namespace Algorithms
continue;
}

// - check peak's resolution
double widthdevpos = width/centre;
if (m_hasInputResolution)
{
double recres = m_inputResolutionWS->readY(wi)[0];
double resmax = recres*m_maxResFactor;
double resmin = recres*m_minResFactor;
if (widthdevpos < resmin || widthdevpos > resmax)
{
std::stringstream dbss;
dbss << " wi = " << wi << " c = " << centre << " Delta(d)/d = " << widthdevpos
<< " too far away from suggested value " << recres;
g_log.debug(dbss.str());
continue;
}
}

// background value
double back_intercept = peakslist->getRef<double>("backgroundintercept", i);
double back_slope = peakslist->getRef<double>("backgroundslope", i);
Expand All @@ -1090,7 +1119,6 @@ namespace Algorithms
if (height < 0.5 * std::sqrt(height + background))
continue;


// - calcualte offsets as to determine the (z-value)
double offset = fabs(peakPositionRef[i]/centre - 1);
if (offset > m_maxOffset)
Expand All @@ -1103,7 +1131,7 @@ namespace Algorithms
else vec_offsets.push_back(offset);

// (g) calculate width/pos as to determine the (z-value) for constant "width" - (delta d)/d
double widthdevpos = width/centre;
// double widthdevpos = width/centre;
vec_widthDivPos.push_back(widthdevpos);

// g_log.debug() << " h:" << height << " c:" << centre << " w:" << (width/(2.*std::sqrt(2.*std::log(2.))))
Expand Down
119 changes: 119 additions & 0 deletions Code/Mantid/Framework/Algorithms/test/GetDetOffsetsMultiPeaksTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,125 @@ class GetDetOffsetsMultiPeaksTest : public CxxTest::TestSuite
TS_ASSERT( !mask->getInstrument()->getDetector(1)->isMasked() );
}

//----------------------------------------------------------------------------------------------
/** Test using the resolution workspace as input
*/
void testExecInputResolutionWS()
{
// ---- Create the simple workspace -------
// Data workspace
MatrixWorkspace_sptr WS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(1,200);
AnalysisDataService::Instance().addOrReplace("temp_event_ws", WS);
WS->getAxis(0)->unit() = Mantid::Kernel::UnitFactory::Instance().create("dSpacing");
const Mantid::MantidVec &X = WS->readX(0);
Mantid::MantidVec &Y = WS->dataY(0);
Mantid::MantidVec &E = WS->dataE(0);
for (size_t i = 0; i < Y.size(); ++i)
{
const double x = (X[i]+X[i+1])/2;
Y[i] = 5.1*exp(-0.5*pow((x-10)/1.0,2));
E[i] = 0.001;
}

// Resolution workspace
MatrixWorkspace_sptr resWS = WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1);
resWS->dataY(0)[0] = 0.2;
AnalysisDataService::Instance().addOrReplace("temp_res_ws", resWS);

// ---- Run algo -----
if ( !offsets.isInitialized() ) offsets.initialize();
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("InputWorkspace","temp_event_ws" ) );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("InputResolutionWorkspace", "temp_res_ws") );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("MinimumResolutionFacyor", 0.8) );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("MinimumResolutionFacyor", 1.2) );

std::string outputWS("offsetsped");
std::string maskWS("masksped");
TS_ASSERT_THROWS_NOTHING( offsets.setPropertyValue("OutputWorkspace",outputWS) );
TS_ASSERT_THROWS_NOTHING( offsets.setPropertyValue("MaskWorkspace",maskWS) );
TS_ASSERT_THROWS_NOTHING(offsets.setPropertyValue("DReference","9.98040"));
TS_ASSERT_THROWS_NOTHING(offsets.setPropertyValue("SpectraFitInfoTableWorkspace", "FitInfoTable"));
TS_ASSERT_THROWS_NOTHING( offsets.execute() );
TS_ASSERT( offsets.isExecuted() );

MatrixWorkspace_const_sptr output;
TS_ASSERT_THROWS_NOTHING( output = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(outputWS) );
if (!output) return;

TS_ASSERT_DELTA( output->dataY(0)[0], -0.002, 0.0002);

AnalysisDataService::Instance().remove(outputWS);

MatrixWorkspace_const_sptr mask;
TS_ASSERT_THROWS_NOTHING( mask = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(maskWS) );
if (!mask) return;
TS_ASSERT( !mask->getInstrument()->getDetector(1)->isMasked() );
}

//----------------------------------------------------------------------------------------------
/** Test using the resolution workspace as input with a failure case
* in which the data is noisy and not valid peak can be found
*/
void testFailInputResolutionWS()
{
// ---- Create the simple workspace -------
// Data workspace
MatrixWorkspace_sptr WS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(1,200);
AnalysisDataService::Instance().addOrReplace("temp_event_ws", WS);
WS->getAxis(0)->unit() = Mantid::Kernel::UnitFactory::Instance().create("dSpacing");
generateNoisyData(WS);

// Resolution workspace
MatrixWorkspace_sptr resWS = WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1);
resWS->dataY(0)[0] = 0.2;
AnalysisDataService::Instance().addOrReplace("temp_res_ws", resWS);

// ---- Run algo -----
if ( !offsets.isInitialized() ) offsets.initialize();
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("InputWorkspace","temp_event_ws" ) );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("InputResolutionWorkspace", "temp_res_ws") );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("MinimumResolutionFacyor", 0.8) );
TS_ASSERT_THROWS_NOTHING( offsets.setProperty("MinimumResolutionFacyor", 1.2) );

std::string outputWS("offsetsped");
std::string maskWS("masksped");
TS_ASSERT_THROWS_NOTHING( offsets.setPropertyValue("OutputWorkspace",outputWS) );
TS_ASSERT_THROWS_NOTHING( offsets.setPropertyValue("MaskWorkspace",maskWS) );
TS_ASSERT_THROWS_NOTHING(offsets.setPropertyValue("DReference","9.98040"));
TS_ASSERT_THROWS_NOTHING(offsets.setPropertyValue("SpectraFitInfoTableWorkspace", "FitInfoTable"));
TS_ASSERT_THROWS_NOTHING( offsets.execute() );
TS_ASSERT( offsets.isExecuted() );

MatrixWorkspace_const_sptr output;
TS_ASSERT_THROWS_NOTHING( output = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(outputWS) );
if (!output) return;

TS_ASSERT_DELTA( output->dataY(0)[0], -0.002, 0.0002);

AnalysisDataService::Instance().remove(outputWS);

MatrixWorkspace_const_sptr mask;
TS_ASSERT_THROWS_NOTHING( mask = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(maskWS) );
if (!mask) return;
TS_ASSERT( !mask->getInstrument()->getDetector(1)->isMasked() );
}

//----------------------------------------------------------------------------------------------
/** Generate noisy data in a workspace
*/
void generateNoisyData(MatrixWorkspace_sptr WS)
{
const Mantid::MantidVec &X = WS->readX(0);
Mantid::MantidVec &Y = WS->dataY(0);
Mantid::MantidVec &E = WS->dataE(0);
for (size_t i = 0; i < Y.size(); ++i)
{
Y[i] = static_cast<double>(rand()%5);
E[i] = 0.01;
}

return;
}

private:
GetDetOffsetsMultiPeaks offsets;
Expand Down

0 comments on commit b968c07

Please sign in to comment.