Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/bugfix/8582_evs_results_no_backg…
Browse files Browse the repository at this point in the history 
…round'
  • Loading branch information
Samuel Jackson committed Dec 10, 2013
2 parents 5fb852c + 065de67 commit 620b652
Show file tree
Hide file tree
Showing 5 changed files with 209 additions and 57 deletions.
10 changes: 7 additions & 3 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@

#include "MantidCurveFitting/ComptonProfile.h"

#include <boost/unordered_map.hpp>

namespace Mantid
{
namespace CurveFitting
Expand Down Expand Up @@ -57,11 +59,11 @@ namespace Mantid
};

/// Calculate & correct the given index of the input workspace
void applyCorrection(const size_t wsIndex);
void applyCorrection(const size_t inputIndex,const size_t outputIndex);
/// Compute the expected spectrum from a given detector
void calculateSpectrumFromDetector(const size_t wsIndex);
void calculateSpectrumFromDetector(const size_t inputIndex, const size_t outputIndex);
/// Compute the expected background from the foils
void calculateBackgroundFromFoils(const size_t wsIndex);
void calculateBackgroundFromFoils(const size_t inputIndex, const size_t outputIndex);
/// Compute expected background from single foil for spectrum at wsIndex
void calculateBackgroundSingleFoil(std::vector<double> & ctfoil,const size_t wsIndex,
const FoilInfo & foilInfo, const Kernel::V3D & detPos,
Expand All @@ -84,6 +86,8 @@ namespace Mantid

/// Input TOF data
API::MatrixWorkspace_const_sptr m_inputWS;
/// Sorted indices to correct
boost::unordered_map<size_t,size_t> m_indices;
/// Function that defines the mass profile
std::string m_profileFunction;
/// The number of peaks in spectrum
Expand Down
105 changes: 71 additions & 34 deletions Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ namespace Mantid

/// Default constructor
CalculateGammaBackground::CalculateGammaBackground()
: Algorithm(), m_inputWS(), m_profileFunction(), m_npeaks(0), m_reversed(),
: Algorithm(), m_inputWS(), m_indices(), m_profileFunction(), m_npeaks(0), m_reversed(),
m_samplePos(), m_l1(0.0), m_foilRadius(0.0), m_foilUpMin(0.0),m_foilUpMax(0.0), m_foils0(), m_foils1(),
m_backgroundWS(), m_correctedWS(), m_progress(NULL)
{
Expand Down Expand Up @@ -102,6 +102,10 @@ namespace Mantid
"Function that is able to compute the mass spectrum for the input data"
"This will usually be the output from the Fitting");

declareProperty(new ArrayProperty<int>("WorkspaceIndexList"),
"Indices of the spectra to include in the correction. If provided, the output only include these spectra\n"
"(Default: all spectra from input)");

declareProperty(new WorkspaceProperty<>("BackgroundWorkspace", "", Direction::Output));
declareProperty(new WorkspaceProperty<>("CorrectedWorkspace", "", Direction::Output));
}
Expand All @@ -111,7 +115,7 @@ namespace Mantid
retrieveInputs();
createOutputWorkspaces();

const int64_t nhist = static_cast<int64_t>(m_inputWS->getNumberHistograms());
const int64_t nhist = static_cast<int64_t>(m_indices.size());
const int64_t nreports = 10 + nhist*(m_npeaks + 2*m_foils0.size()*NTHETA*NUP*m_npeaks);
m_progress = new Progress(this, 0.0, 1.0, nreports);

Expand All @@ -120,25 +124,34 @@ namespace Mantid
{
PARALLEL_START_INTERUPT_REGION

m_backgroundWS->setX(i,m_inputWS->refX(i));
m_correctedWS->setX(i,m_inputWS->refX(i));
const size_t outputIndex = i;
auto indexIter = m_indices.cbegin();
std::advance(indexIter, i);
const size_t inputIndex = indexIter->second;

m_backgroundWS->setX(outputIndex,m_inputWS->refX(inputIndex));
m_correctedWS->setX(outputIndex,m_inputWS->refX(inputIndex));
try
{
specid_t spectrumNo = m_inputWS->getSpectrum(i)->getSpectrumNo();
const auto * inSpec = m_inputWS->getSpectrum(inputIndex);
const specid_t spectrumNo(inSpec->getSpectrumNo());
m_backgroundWS->getSpectrum(outputIndex)->copyInfoFrom(*inSpec);
m_correctedWS->getSpectrum(outputIndex)->copyInfoFrom(*inSpec);

if(spectrumNo >= FORWARD_SCATTER_SPECMIN && spectrumNo <= FORWARD_SCATTER_SPECMAX )
{
applyCorrection(i);
applyCorrection(inputIndex,outputIndex);
}
else
{
g_log.information("Spectrum " + boost::lexical_cast<std::string>(spectrumNo) + " not in forward scatter range. Skipping correction.");
// Leave background at 0 and just copy data to corrected
m_correctedWS->dataY(i) = m_inputWS->readY(i);
m_correctedWS->dataY(outputIndex) = m_inputWS->readY(inputIndex);
}
}
catch(Exception::NotFoundError &)
{
g_log.information("No detector defined for index=" + boost::lexical_cast<std::string>(i) + ". Skipping correction.");
g_log.information("No detector defined for index=" + boost::lexical_cast<std::string>(inputIndex) + ". Skipping correction.");
}

PARALLEL_END_INTERUPT_REGION
Expand All @@ -153,26 +166,28 @@ namespace Mantid
/**
* Calculate & apply gamma correction for the given index of the
* input workspace
* @param wsIndex A workspace index that defines the spectrum to correct
* @param inputIndex A workspace index that defines the input spectrum to correct
* @param outputIndex A workspace index that defines the output to hold the results
*/
void CalculateGammaBackground::applyCorrection(const size_t wsIndex)
void CalculateGammaBackground::applyCorrection(const size_t inputIndex, const size_t outputIndex)
{
m_progress->report("Computing TOF from detector");

// results go straight in m_correctedWS to save memory allocations
calculateSpectrumFromDetector(wsIndex);
calculateSpectrumFromDetector(inputIndex,outputIndex);

m_progress->report("Computing TOF foils");
// Output goes to m_background to save memory allocations
calculateBackgroundFromFoils(wsIndex);
calculateBackgroundFromFoils(inputIndex,outputIndex);

m_progress->report("Computing correction to input");
// Compute total counts from input data, (detector-foil) contributions
//assume constant binning
const size_t nbins = m_correctedWS->blocksize();
const auto & inY = m_inputWS->readY(wsIndex);
auto & detY = m_correctedWS->dataY(wsIndex);
auto & foilY = m_backgroundWS->dataY(wsIndex);
const double deltaT = m_correctedWS->readX(wsIndex)[1] - m_correctedWS->readX(wsIndex)[0];
const auto & inY = m_inputWS->readY(inputIndex);
auto & detY = m_correctedWS->dataY(outputIndex);
auto & foilY = m_backgroundWS->dataY(outputIndex);
const double deltaT = m_correctedWS->readX(outputIndex)[1] - m_correctedWS->readX(outputIndex)[0];

double dataCounts(0.0), simulCounts(0.0);
for(size_t j = 0; j < nbins; ++j)
Expand All @@ -196,12 +211,13 @@ namespace Mantid
}

/**
* Results are placed in the corresponding index on the output corrected workspace
* @param wsIndex The spectrum to compute
* Results are placed in the mapped index on the output corrected workspace
* @param inputIndex Workspace index that defines the input spectrum to correct
* @param outputIndex Workspace index that defines the spectrum to hold the results
*/
void CalculateGammaBackground::calculateSpectrumFromDetector(const size_t wsIndex)
void CalculateGammaBackground::calculateSpectrumFromDetector(const size_t inputIndex, const size_t outputIndex)
{
auto det = m_inputWS->getDetector(wsIndex);
auto det = m_inputWS->getDetector(inputIndex);
auto detPos = det->getPos();

// -- Setup detector & resolution parameters --
Expand All @@ -221,9 +237,9 @@ namespace Mantid

// Compute a time of flight spectrum convolved with a Voigt resolution function for each mass
// at the detector point & sum to a single spectrum
auto & ctdet = m_correctedWS->dataY(wsIndex);
auto & ctdet = m_correctedWS->dataY(outputIndex);
std::vector<double> tmpWork(ctdet.size());
calculateTofSpectrum(ctdet,tmpWork, wsIndex, detPar, detRes);
calculateTofSpectrum(ctdet,tmpWork, outputIndex, detPar, detRes);
//Correct for distance to the detector: 0.5/l2^2
const double detDistCorr = 0.5/detPar.l2/detPar.l2;
std::transform(ctdet.begin(),ctdet.end(),ctdet.begin(),
Expand All @@ -233,11 +249,12 @@ namespace Mantid
/**
* Calculate & apply gamma correction for the given index of the
* input workspace
* @param wsIndex A workspace index that defines the spectrum to correct
* @param inputIndex Workspace index that defines the input spectrum to correct
* @param outputIndex Workspace index that defines the spectrum to hold the results
*/
void CalculateGammaBackground::calculateBackgroundFromFoils(const size_t wsIndex)
void CalculateGammaBackground::calculateBackgroundFromFoils(const size_t inputIndex, const size_t outputIndex)
{
auto det = m_inputWS->getDetector(wsIndex);
auto det = m_inputWS->getDetector(inputIndex);
auto detPos = det->getPos();

// -- Setup detector & resolution parameters --
Expand All @@ -257,25 +274,25 @@ namespace Mantid

const size_t nxvalues = m_backgroundWS->blocksize();
std::vector<double> foilSpectrum(nxvalues);
auto & ctfoil = m_backgroundWS->dataY(wsIndex);
auto & ctfoil = m_backgroundWS->dataY(outputIndex);

// Compute (C1 - C0) where C1 is counts in pos 1 and C0 counts in pos 0
assert(m_foils0.size() == m_foils1.size());
for(size_t i = 0; i < m_foils0.size(); ++i)
{
foilSpectrum.assign(nxvalues,0.0);
calculateBackgroundSingleFoil(foilSpectrum, wsIndex,m_foils1[i], detPos, detPar, detRes);
calculateBackgroundSingleFoil(foilSpectrum, outputIndex,m_foils1[i], detPos, detPar, detRes);
// sum spectrum values from first position
std::transform(ctfoil.begin(), ctfoil.end(), foilSpectrum.begin(), ctfoil.begin(),
std::plus<double>());

foilSpectrum.assign(nxvalues,0.0);
calculateBackgroundSingleFoil(foilSpectrum, wsIndex,m_foils0[i], detPos, detPar, detRes);
calculateBackgroundSingleFoil(foilSpectrum, outputIndex,m_foils0[i], detPos, detPar, detRes);
// subtract spectrum values from zeroth position
std::transform(ctfoil.begin(), ctfoil.end(), foilSpectrum.begin(), ctfoil.begin(),
std::minus<double>());
}
bool reversed = (m_reversed.count(m_inputWS->getSpectrum(wsIndex)->getSpectrumNo()));
bool reversed = (m_reversed.count(m_inputWS->getSpectrum(inputIndex)->getSpectrumNo()));
// This is quicker than the if within the loop
if(reversed)
{
Expand All @@ -290,7 +307,7 @@ namespace Mantid
* Integrates over the foil area defined by the foil radius to accumulate an estimate of the counts
* resulting from this region
* @param ctfoil Output vector to hold results
* @param wsIndex Index on workspaces currently operating on
* @param wsIndex Index on output background workspaces currently operating
* @param foilInfo Foil description object
* @param detPos The pre-calculated detector V3D
* @param detPar DetectorParams object that defines information on the detector associated with spectrum at wsIndex
Expand Down Expand Up @@ -356,7 +373,7 @@ namespace Mantid
* Uses the compton profile functions to compute a particular mass spectrum
* @param result [Out] The value of the computed spectrum
* @param tmpWork [In] Pre-allocated working area that will be overwritten
* @param wsIndex Index on the input workspace that defines the detector parameters
* @param wsIndex Index on the output background workspace that gives the X values to use
* @param detpar Struct containing parameters about the detector
* @param respar Struct containing parameters about the resolution
*/
Expand Down Expand Up @@ -424,7 +441,6 @@ namespace Mantid
"composite of ComptonProfiles or a single ComptonProfile.");
}

cacheInstrumentGeometry();
// Spectrum numbers whose calculation of background from foils is reversed
m_reversed.clear();
for(specid_t i = 143; i < 199; ++i)
Expand All @@ -433,6 +449,26 @@ namespace Mantid
(i >= 175 && i <= 182) || (i >= 191 && i <= 198) )
m_reversed.insert(i);
}

// Workspace indices mapping input->output
m_indices.clear();
std::vector<int> requestedIndices = getProperty("WorkspaceIndexList");
if(requestedIndices.empty())
{
for(size_t i = 0; i < m_inputWS->getNumberHistograms(); ++i)
{
m_indices[i] = i; // 1-to-1
}
}
else
{
for(size_t i = 0; i < requestedIndices.size(); ++i)
{
m_indices[i] = static_cast<size_t>(requestedIndices[i]); //user-requested->increasing on output
}
}

cacheInstrumentGeometry();
}


Expand All @@ -441,8 +477,9 @@ namespace Mantid
*/
void CalculateGammaBackground::createOutputWorkspaces()
{
m_backgroundWS = WorkspaceFactory::Instance().create(m_inputWS);
m_correctedWS = WorkspaceFactory::Instance().create(m_inputWS);
const size_t nhist = m_indices.size();
m_backgroundWS = WorkspaceFactory::Instance().create(m_inputWS,nhist);
m_correctedWS = WorkspaceFactory::Instance().create(m_inputWS,nhist);
}

/**
Expand Down
71 changes: 59 additions & 12 deletions Code/Mantid/Framework/CurveFitting/src/ComptonScatteringCountRate.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "MantidCurveFitting/ComptonScatteringCountRate.h"
#include "MantidCurveFitting/AugmentedLagrangianOptimizer.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidKernel/Math/Optimization/SLSQPMinimizer.h"

#include <boost/bind.hpp>

Expand Down Expand Up @@ -64,11 +65,42 @@ namespace CurveFitting

namespace
{
/// Defines a helper struct to compute 0.5*||Cx-d||^2
struct Norm2
/// Struct to compute norm2,0.5*||Cx-d||^2, when the background is NOT included and the
/// optimizer expects constraints specified as >= 0
struct NoBkgdNorm2
{
/// Compute the value of the objective function
Norm2(const Kernel::DblMatrix & cmatrix, const std::vector<double> & data)
NoBkgdNorm2(const Kernel::DblMatrix & cmatrix, const std::vector<double> & data)
: cm(cmatrix), nrows(cmatrix.numRows()), ncols(cmatrix.numCols()), rhs(data) {}

double eval(const std::vector<double> & xpt) const
{
double norm2(0.0);
for(size_t i = 0; i < nrows; ++i)
{
const double *cmRow = cm[i];
double cx(0.0);
for(size_t j = 0; j < ncols; ++j)
{
cx += cmRow[j]*xpt[j];
}
cx -= rhs[i];
norm2 += cx*cx;
}
return 0.5*norm2;
}
const Kernel::DblMatrix & cm;
size_t nrows;
size_t ncols;
const std::vector<double> & rhs;
};
/// Struct to compute norm2 when the background is included and the
/// optimizer expects constraints specified as <= 0.
/// It is assumed the CM matrix has been multiplied by -1
struct BkgdNorm2
{
/// Compute the value of the objective function
BkgdNorm2(const Kernel::DblMatrix & cmatrix, const std::vector<double> & data)
: cm(cmatrix), nrows(cmatrix.numRows()), ncols(cmatrix.numCols()), rhs(data) {}

double eval(const size_t n, const double * xpt) const
Expand All @@ -93,6 +125,7 @@ namespace CurveFitting
size_t ncols;
const std::vector<double> & rhs;
};

}

/**
Expand Down Expand Up @@ -120,14 +153,23 @@ namespace CurveFitting
// Compute the constraint matrix
this->updateCMatrixValues();

Norm2 objf(m_cmatrix, m_dataErrorRatio);
//boost::function<double(const size_t,const double *)> objfunc =
AugmentedLagrangianOptimizer::ObjFunction objfunc = boost::bind(&Norm2::eval, objf, _1, _2);
AugmentedLagrangianOptimizer lsqmin(nparams, objfunc, m_eqMatrix, m_cmatrix);
lsqmin.minimize(x0);

// Set the parameters for the 'real' function calls
setFixedParameterValues(x0);
if(m_bkgdPolyN > 0)
{
BkgdNorm2 objf(m_cmatrix, m_dataErrorRatio);
AugmentedLagrangianOptimizer::ObjFunction objfunc = boost::bind(&BkgdNorm2::eval, objf, _1, _2);
AugmentedLagrangianOptimizer lsqmin(nparams, objfunc, m_eqMatrix, m_cmatrix);
lsqmin.minimize(x0);
// Set the parameters for the 'real' function calls
setFixedParameterValues(x0);
}
else
{
NoBkgdNorm2 objfunc(m_cmatrix, m_dataErrorRatio);
Kernel::Math::SLSQPMinimizer lsqmin(nparams, objfunc, m_eqMatrix, m_cmatrix);
auto res = lsqmin.minimize(x0);
// Set the parameters for the 'real' function calls
setFixedParameterValues(res);
}
}

/**
Expand Down Expand Up @@ -167,7 +209,12 @@ namespace CurveFitting
const size_t numFilled = profile->fillConstraintMatrix(m_cmatrix,start,m_errors);
start += numFilled;
}
m_cmatrix *= -1.0;
// Using different minimizer for background constraints as the original is not stable
// The background one requires the contraints are specified as <= 0
if(m_bkgdPolyN > 0)
{
m_cmatrix *= -1.0;
}

if(g_log.is(Logger::Priority::PRIO_DEBUG))
{
Expand Down

0 comments on commit 620b652

Please sign in to comment.