Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Refactored FindPeakBackground. Refs #9020.
1. Refactored and thus fixed some issues; 2. Added more unit tests; 3. Change input property. Algorithm can work with 1 spectrum only.
- Loading branch information
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -76,9 +76,8 @@ namespace Algorithms | |
| auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input); | ||
| declareProperty(inwsprop, "Name of input MatrixWorkspace that contains peaks."); | ||
|
|
||
| declareProperty("WorkspaceIndex", EMPTY_INT(), "workspace indices to have peak and background separated. " | ||
| "No default is taken. "); | ||
|
|
||
| declareProperty("SigmaConstant", 1.0, "Multiplier of standard deviations of the variance for convergence of " | ||
| "peak elimination. Default is 1.0. "); | ||
|
|
@@ -109,32 +108,27 @@ namespace Algorithms | |
| { | ||
| // 1. Get input and validate | ||
| MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace"); | ||
| int inpwsindex = getProperty("WorkspaceIndex"); | ||
| std::vector<double> m_vecFitWindows = getProperty("FitWindow"); | ||
| m_backgroundType = getPropertyValue("BackgroundType"); | ||
| double k = getProperty("SigmaConstant"); | ||
|
|
||
| if (inpwsindex < 0 || inpwsindex >= static_cast<int>(inpWS->getNumberHistograms())) | ||
| { | ||
| stringstream errss; | ||
| errss << "Input workspace " << inpWS->name() << " has " << inpWS->getNumberHistograms() | ||
| << " spectra. Input workspace index " << inpwsindex << " is out of boundary. "; | ||
| throw runtime_error(errss.str()); | ||
| } | ||
| if (isEmpty(inpwsindex)) | ||
| throw runtime_error("WorkspaceIndex must be given. "); | ||
|
|
||
| // 2. Generate output | ||
| const MantidVec& inpX = inpWS->readX(inpwsindex); | ||
| size_t sizex = inpWS->readX(inpwsindex).size(); | ||
| size_t sizey = inpWS->readY(inpwsindex).size(); | ||
| size_t n = sizey; | ||
| size_t l0 = 0; | ||
|
|
||
| if (m_vecFitWindows.size() > 1) | ||
| { | ||
|
|
@@ -152,145 +146,121 @@ namespace Algorithms | |
| m_outPeakTableWS->addColumn("double", "bkg0"); | ||
| m_outPeakTableWS->addColumn("double", "bkg1"); | ||
| m_outPeakTableWS->addColumn("double", "bkg2"); | ||
|
|
||
| m_outPeakTableWS->appendRow(); | ||
|
|
||
| // 3. Get Y values | ||
| Progress prog(this, 0, 1.0, 1); | ||
This comment has been minimized.
Sorry, something went wrong.
This comment has been minimized.
Sorry, something went wrong. 
RussellTaylor
via email
Contributor
|
||
|
|
||
| // Find background | ||
|
|
||
| const MantidVec& inpY = inpWS->readY(inpwsindex); | ||
|
|
||
| double Ymean, Yvariance, Ysigma; | ||
| MantidVec maskedY; | ||
| MantidVec::const_iterator in = std::min_element(inpY.begin(), inpY.end()); | ||
| double bkg0 = inpY[in - inpY.begin()]; | ||
| for (size_t l = l0; l < n; ++l) | ||
| { | ||
| maskedY.push_back(inpY[l]-bkg0); | ||
| } | ||
| MantidVec mask(n-l0,0.0); | ||
| double xn = static_cast<double>(n-l0); | ||
| do | ||
| { | ||
| Statistics stats = getStatistics(maskedY); | ||
| Ymean = stats.mean; | ||
| Yvariance = stats.standard_deviation * stats.standard_deviation; | ||
| Ysigma = std::sqrt((moment4(maskedY,n-l0,Ymean)-(xn-3.0)/(xn-1.0) * Yvariance)/xn); | ||
| MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end()); | ||
| const size_t pos = it - maskedY.begin(); | ||
| maskedY[pos] = 0; | ||
| mask[pos] = 1.0; | ||
| } | ||
| while (std::abs(Ymean-Yvariance) > k * Ysigma); | ||
|
|
||
| if(n-l0 > 5) | ||
| { | ||
| // remove single outliers | ||
| if (mask[1] == mask[2] && mask[2] == mask[3]) | ||
| mask[0] = mask[1]; | ||
| if (mask[0] == mask[2] && mask[2] == mask[3]) | ||
| mask[1] = mask[2]; | ||
| for (size_t l = 2; l < n-l0-3; ++l) | ||
| { | ||
| if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2])) | ||
| { | ||
| mask[l] = mask[l+1]; | ||
| } | ||
| } | ||
| if (mask[n-l0-2] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4]) | ||
| mask[n-l0-1] = mask[n-l0-2]; | ||
| if (mask[n-l0-1] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4]) | ||
| mask[n-l0-2] = mask[n-l0-1]; | ||
|
|
||
| // mask regions not connected to largest region | ||
| // for loop can start > 1 for multiple peaks | ||
| vector<cont_peak> peaks; | ||
| if (mask[0] == 1) | ||
| { | ||
| peaks.push_back(cont_peak()); | ||
| peaks[peaks.size()-1].start = l0; | ||
| } | ||
| for (size_t l = 1; l < n-l0; ++l) | ||
| { | ||
| if (mask[l] != mask[l-1] && mask[l] == 1) | ||
| { | ||
| peaks.push_back(cont_peak()); | ||
| peaks[peaks.size()-1].start = l+l0; | ||
| } | ||
| else if (peaks.size() > 0) | ||
| { | ||
| size_t ipeak = peaks.size()-1; | ||
| if (mask[l] != mask[l-1] && mask[l] == 0) | ||
| { | ||
| peaks[ipeak].stop = l+l0; | ||
| } | ||
| if (inpY[l+l0] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l+l0]; | ||
| } | ||
| } | ||
| size_t min_peak, max_peak; | ||
| double a0,a1,a2; | ||
| if(peaks.size()> 0) | ||
| { | ||
| if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1; | ||
| std::sort(peaks.begin(), peaks.end(), by_len()); | ||
|
|
||
| // save endpoints | ||
| min_peak = peaks[0].start; | ||
| // extra point for histogram input | ||
| max_peak = peaks[0].stop + sizex - sizey; | ||
| estimateBackground(inpX, inpY, l0, n, | ||
| peaks[0].start, peaks[0].stop, a0, a1, a2); | ||
| } | ||
| else | ||
| { | ||
| // assume background is 12 first and last points | ||
| min_peak = l0+12; | ||
| max_peak = n-13; | ||
| if (min_peak > sizey)min_peak = sizey-1; | ||
| a0 = 0.0; | ||
| a1 = 0.0; | ||
| a2 = 0.0; | ||
| } | ||
|
|
||
| // Add a new row | ||
| API::TableRow t = m_outPeakTableWS->getRow(0); | ||
| t << static_cast<int>(inpwsindex) << static_cast<int>(min_peak) << static_cast<int>(max_peak) << a0 << a1 <<a2; | ||
| } | ||
|
|
||
| prog.report(); | ||
|
|
||
| // 4. Set the output | ||
| setProperty("OutputWorkspace", m_outPeakTableWS); | ||
|
|
||
| return; | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------------------------- | ||
| /** Estimate background | ||
| * @param X :: vec for X | ||
|
|
||
Oops, something went wrong.
This line is what is causing the error on the autoreduction. There is a check in Algorithm.cpp (around line 817) that checks
if (startProgress >= 0 && endProgress > startProgress && endProgress <= 1.). When that check is true, an error is written out to the log, which in turn makes the autoreduction report an error. That is also why the reduction completed, since that check doesn't actually raise any exception.