Unverified Commit faac501e authored by Savici, Andrei T's avatar Savici, Andrei T Committed by GitHub
Browse files

Merge pull request #23674 from mantidproject/23614_peak_integration_improvements

Peak integration improvements
parents 718000e6 d4b5267d
......@@ -633,7 +633,7 @@ TMDE(void MDBox)::integrateCylinder(
for (const auto &evnt : events) {
coord_t out[2]; // radius and length of cylinder
radiusTransform.apply(evnt.getCenter(), out);
if (out[0] < radius && std::fabs(out[1]) < 0.5 * length) {
if (out[0] < radius && std::fabs(out[1]) < 0.5 * length + deltaQ) {
// add event to appropriate y channel
size_t xchannel =
static_cast<size_t>(std::floor(out[1] / deltaQ)) + numSteps / 2;
......
......@@ -489,71 +489,56 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
}
errorSquared = std::fabs(signal);
} else {
API::IAlgorithm_sptr findpeaks =
createChildAlgorithm("FindPeaks", -1, -1, false);
findpeaks->setProperty("InputWorkspace", wsProfile2D);
findpeaks->setProperty<int>("FWHM", 7);
findpeaks->setProperty<int>("Tolerance", 4);
// FindPeaks will do the checking on the validity of WorkspaceIndex
findpeaks->setProperty("WorkspaceIndex", static_cast<int>(i));
// Get the specified peak positions, which is optional
findpeaks->setProperty<std::string>("PeakFunction", profileFunction);
// FindPeaks will use linear or flat if they are better
findpeaks->setProperty<std::string>("BackgroundType", "Quadratic");
findpeaks->setProperty<bool>("HighBackground", true);
findpeaks->setProperty<bool>("RawPeakParameters", true);
std::vector<double> peakPosToFit;
peakPosToFit.push_back(static_cast<double>(numSteps) / 2.0);
findpeaks->setProperty("PeakPositions", peakPosToFit);
findpeaks->setProperty<int>("MinGuessedPeakWidth", 4);
findpeaks->setProperty<int>("MaxGuessedPeakWidth", 4);
IAlgorithm_sptr fitAlgorithm =
createChildAlgorithm("Fit", -1, -1, false);
// fitAlgorithm->setProperty("CreateOutput", true);
// fitAlgorithm->setProperty("Output", "FitPeaks1D");
std::string myFunc =
std::string("name=LinearBackground;name=") + profileFunction;
auto maxPeak = std::max_element(signal_fit.begin(), signal_fit.end());
std::ostringstream strs;
strs << maxPeak[0];
std::string strMax = strs.str();
if (profileFunction == "Gaussian") {
myFunc += ", PeakCentre=50, Height=" + strMax;
fitAlgorithm->setProperty("Constraints", "40<f1.PeakCentre<60");
} else if (profileFunction == "BackToBackExponential" ||
profileFunction == "IkedaCarpenterPV") {
myFunc += ", X0=50, I=" + strMax;
fitAlgorithm->setProperty("Constraints", "40<f1.X0<60");
}
fitAlgorithm->setProperty("CalcErrors", true);
fitAlgorithm->setProperty("Function", myFunc);
fitAlgorithm->setProperty("InputWorkspace", wsProfile2D);
fitAlgorithm->setProperty("WorkspaceIndex", static_cast<int>(i));
try {
findpeaks->executeAsChildAlg();
fitAlgorithm->executeAsChildAlg();
} catch (...) {
g_log.error("Can't execute FindPeaks algorithm");
g_log.error("Can't execute Fit algorithm");
continue;
}
API::ITableWorkspace_sptr paramws = findpeaks->getProperty("PeaksList");
if (paramws->rowCount() < 1)
continue;
std::ostringstream fun_str;
fun_str << "name=" << profileFunction;
size_t numcols = paramws->columnCount();
std::vector<std::string> paramsName = paramws->getColumnNames();
std::vector<double> paramsValue;
API::TableRow row = paramws->getRow(0);
int spectrum;
row >> spectrum;
for (size_t j = 1; j < numcols; ++j) {
double parvalue;
row >> parvalue;
if (j == numcols - 4)
fun_str << ";name=Quadratic";
// erase f0. or f1.
// if (j > 0 && j < numcols-1) fun_str << "," <<
// paramsName[j].erase(0,3) <<"="<<parvalue;
if (j > 0 && j < numcols - 1)
fun_str << "," << paramsName[j] << "=" << parvalue;
paramsValue.push_back(parvalue);
}
IFunction_sptr ifun = fitAlgorithm->getProperty("Function");
if (i == 0) {
for (size_t j = 0; j < numcols; ++j)
out << std::setw(20) << paramsName[j] << " ";
out << std::setw(20) << "spectrum"
<< " ";
for (size_t j = 0; j < ifun->nParams(); ++j)
out << std::setw(20) << ifun->parameterName(j) << " ";
out << std::setw(20) << "chi2"
<< " ";
out << "\n";
}
out << std::setw(20) << i;
for (size_t j = 0; j < numcols - 1; ++j)
out << std::setw(20) << i << " ";
for (size_t j = 0; j < ifun->nParams(); ++j) {
out << std::setw(20) << std::fixed << std::setprecision(10)
<< paramsValue[j] << " ";
out << "\n";
// Evaluate fit at points
<< ifun->getParameter(j) << " ";
}
double chi2 = fitAlgorithm->getProperty("OutputChi2overDoF");
out << std::setw(20) << std::fixed << std::setprecision(10) << chi2
<< "\n";
IFunction_sptr ifun =
FunctionFactory::Instance().createInitialized(fun_str.str());
boost::shared_ptr<const CompositeFunction> fun =
boost::dynamic_pointer_cast<const CompositeFunction>(ifun);
......@@ -593,10 +578,8 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
errorSquared = std::fabs(signal);
// Get background counts
for (size_t j = 0; j < numSteps; j++) {
// paramsValue[numcols-2] is chisq
double background = paramsValue[numcols - 3] * x[j] * x[j] +
paramsValue[numcols - 4] * x[j] +
paramsValue[numcols - 5];
double background =
ifun->getParameter(0) + ifun->getParameter(1) * x[j];
if (j < peakMin || j > peakMax)
background_total = background_total + background;
}
......
......@@ -33,6 +33,7 @@ Improvements
- :ref:`MDNormSCD <algm-MDNormSCD>` now can handle merged MD workspaces.
- :ref:`StartLiveData <algm-StartLiveData>` will load "live"
data streaming from TOPAZ new Adara data server.
- :ref:`IntegratePeaksMD <algm-IntegratePeaksMD>` with Cylinder=True now has improved fits using BackToBackExponential and IkedaCarpenterPV functions.
Bugfixes
########
......
......@@ -3008,7 +3008,7 @@
</item>
<item>
<property name="text">
<string>Bk2BKExpConvPV</string>
<string>Bk2BkExpConvPV</string>
</property>
</item>
<item>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment