Skip to content

Commit

Permalink
Add support for DetectorGroups to Compton functions
Browse files Browse the repository at this point in the history 
The parameters are averaged across the DetectorGroup when one is found
Refs #8593
  • Loading branch information
@martyngigg
martyngigg committed Dec 10, 2013
1 parent 99e0019 commit 1745c65
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 71 deletions.
2 changes: 0 additions & 2 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h
View file
Expand Up @@ -81,8 +81,6 @@ namespace Mantid
/// Compute the theta range for a given foil
std::pair<double,double> calculateThetaRange(const Geometry::IComponent_const_sptr & foilComp,
const double radius, const unsigned int horizDir) const;
/// Retrieve parameter for given component
double getComponentParameter(const Geometry::IComponent & comp,const std::string &name) const;

/// Input TOF data
API::MatrixWorkspace_const_sptr m_inputWS;
Expand Down
8 changes: 4 additions & 4 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonProfile.h
View file
Expand Up @@ -58,6 +58,10 @@ namespace CurveFitting
class MANTID_CURVEFITTING_DLL ComptonProfile : public virtual API::ParamFunction, public virtual API::IFunction1D
{
public:
/// Retrieve a component parameter
static double getComponentParameter(const Geometry::IComponent_const_sptr & comp, const Geometry::ParameterMap &pmap,
const std::string &name);

/// Default constructor required for factory
ComptonProfile();

Expand Down Expand Up @@ -116,10 +120,6 @@ namespace CurveFitting
/// Logger
Kernel::Logger & m_log;

private:
/// Retrieve a component parameter
double getComponentParameter(const Geometry::IComponent & comp,const std::string &name) const;

/// Current workspace index, required to access instrument parameters
size_t m_wsIndex;
/// Store the mass values
Expand Down
58 changes: 21 additions & 37 deletions Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp
View file
Expand Up @@ -218,22 +218,23 @@ namespace Mantid
void CalculateGammaBackground::calculateSpectrumFromDetector(const size_t inputIndex, const size_t outputIndex)
{
auto det = m_inputWS->getDetector(inputIndex);
const auto & pmap = m_inputWS->constInstrumentParameters();
auto detPos = det->getPos();

// -- Setup detector & resolution parameters --
DetectorParams detPar;
detPar.l1 = m_l1;
detPar.l2 = m_samplePos.distance(detPos);
detPar.theta = m_inputWS->detectorTwoTheta(det); //radians
detPar.t0 = getComponentParameter(*det, "t0")*1e-6; // seconds
detPar.efixed = getComponentParameter(*det,"efixed");
detPar.t0 = ComptonProfile::getComponentParameter(det, pmap, "t0")*1e-6; // seconds
detPar.efixed = ComptonProfile::getComponentParameter(det, pmap,"efixed");

ResolutionParams detRes;
detRes.dl1 = getComponentParameter(*det, "sigma_l1"); // DL0
detRes.dl2 = getComponentParameter(*det, "sigma_l2"); // DL1
detRes.dthe = getComponentParameter(*det, "sigma_theta"); //DTH in radians
detRes.dEnGauss = getComponentParameter(*det, "sigma_gauss");
detRes.dEnLorentz = getComponentParameter(*det, "hwhm_lorentz");
detRes.dl1 = ComptonProfile::getComponentParameter(det, pmap, "sigma_l1"); // DL0
detRes.dl2 = ComptonProfile::getComponentParameter(det, pmap, "sigma_l2"); // DL1
detRes.dthe = ComptonProfile::getComponentParameter(det, pmap, "sigma_theta"); //DTH in radians
detRes.dEnGauss = ComptonProfile::getComponentParameter(det, pmap, "sigma_gauss");
detRes.dEnLorentz = ComptonProfile::getComponentParameter(det, pmap, "hwhm_lorentz");

// Compute a time of flight spectrum convolved with a Voigt resolution function for each mass
// at the detector point & sum to a single spectrum
Expand All @@ -255,22 +256,23 @@ namespace Mantid
void CalculateGammaBackground::calculateBackgroundFromFoils(const size_t inputIndex, const size_t outputIndex)
{
auto det = m_inputWS->getDetector(inputIndex);
const auto & pmap = m_inputWS->constInstrumentParameters();
auto detPos = det->getPos();

// -- Setup detector & resolution parameters --
DetectorParams detPar;
detPar.l1 = m_l1;
detPar.l2 = m_samplePos.distance(detPos);
detPar.theta = m_inputWS->detectorTwoTheta(det); //radians
detPar.t0 = getComponentParameter(*det, "t0")*1e-6; // seconds
detPar.efixed = getComponentParameter(*det,"efixed");
detPar.t0 = ComptonProfile::getComponentParameter(det, pmap, "t0")*1e-6; // seconds
detPar.efixed = ComptonProfile::getComponentParameter(det, pmap,"efixed");

ResolutionParams detRes;
detRes.dl1 = getComponentParameter(*det, "sigma_l1"); // DL0
detRes.dl2 = getComponentParameter(*det, "sigma_l2"); // DL1
detRes.dthe = getComponentParameter(*det, "sigma_theta"); //DTH in radians
detRes.dEnGauss = getComponentParameter(*det, "sigma_gauss");
detRes.dEnLorentz = getComponentParameter(*det, "hwhm_lorentz");
detRes.dl1 = ComptonProfile::getComponentParameter(det, pmap, "sigma_l1"); // DL0
detRes.dl2 = ComptonProfile::getComponentParameter(det, pmap, "sigma_l2"); // DL1
detRes.dthe = ComptonProfile::getComponentParameter(det, pmap, "sigma_theta"); //DTH in radians
detRes.dEnGauss = ComptonProfile::getComponentParameter(det, pmap, "sigma_gauss");
detRes.dEnLorentz = ComptonProfile::getComponentParameter(det, pmap, "hwhm_lorentz");

const size_t nxvalues = m_backgroundWS->blocksize();
std::vector<double> foilSpectrum(nxvalues);
Expand Down Expand Up @@ -510,6 +512,7 @@ namespace Mantid

// foil geometry
// there should be the same number in each position
const auto & pmap = m_inputWS->constInstrumentParameters();
auto foils0 = inst->getAllComponentsWithName("foil-pos0");
auto foils1 = inst->getAllComponentsWithName("foil-pos1");
const size_t nfoils = foils0.size();
Expand All @@ -534,16 +537,16 @@ namespace Mantid
FoilInfo descr;
descr.thetaMin = thetaRng0.first;
descr.thetaMax = thetaRng0.second;
descr.lorentzWidth = getComponentParameter(*foil0, "hwhm_lorentz");
descr.gaussWidth = getComponentParameter(*foil0, "sigma_gauss");
descr.lorentzWidth = ComptonProfile::getComponentParameter(foil0, pmap, "hwhm_lorentz");
descr.gaussWidth = ComptonProfile::getComponentParameter(foil0, pmap, "sigma_gauss");
m_foils0[i] = descr; //copy

const auto & foil1 = foils1[i];
auto thetaRng1 = calculateThetaRange(foil1, m_foilRadius,refFrame->pointingHorizontal());
descr.thetaMin = thetaRng1.first;
descr.thetaMax = thetaRng1.second;
descr.lorentzWidth = getComponentParameter(*foil1, "hwhm_lorentz");
descr.gaussWidth = getComponentParameter(*foil1, "sigma_gauss");
descr.lorentzWidth = ComptonProfile::getComponentParameter(foil1, pmap, "hwhm_lorentz");
descr.gaussWidth = ComptonProfile::getComponentParameter(foil1, pmap, "sigma_gauss");
m_foils1[i] = descr; //copy
}

Expand Down Expand Up @@ -596,24 +599,5 @@ namespace Mantid
return std::make_pair(theta - dtheta, theta + dtheta);
}

/**
* @param comp A reference to the component that should contain the parameter
* @param name The name of the parameter
* @returns The value of the parameter if it exists
* @throws A std::invalid_argument error if the parameter does not exist
*/
double CalculateGammaBackground::getComponentParameter(const Geometry::IComponent & comp,const std::string &name) const
{
std::vector<double> pars = comp.getNumberParameter(name);
if(!pars.empty())
{
return pars[0];
}
else
{
throw std::invalid_argument("CalculateGammaBackground - Cannot find component parameter \"" + name + "\".");
}
}

} // namespace CurveFitting
} // namespace Mantid
78 changes: 50 additions & 28 deletions Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp
View file
Expand Up @@ -3,7 +3,7 @@
//-----------------------------------------------------------------------------
#include "MantidCurveFitting/ComptonProfile.h"
#include "MantidAPI/FunctionFactory.h"

#include "MantidGeometry/Instrument/DetectorGroup.h"
#include <gsl/gsl_poly.h>

namespace Mantid
Expand All @@ -21,6 +21,47 @@ namespace CurveFitting
///@endcond
}

/**
* If a DetectorGroup is encountered then the parameters are averaged over the group
* @param comp A pointer to the component that should contain the parameter
* @param pmap A reference to the ParameterMap that stores the parameters
* @param name The name of the parameter
* @returns The value of the parameter if it exists
* @throws A std::invalid_argument error if the parameter does not exist
*/
double ComptonProfile::getComponentParameter(const Geometry::IComponent_const_sptr & comp, const Geometry::ParameterMap &pmap,
const std::string &name)
{
if(!comp) throw std::invalid_argument("ComptonProfile - Cannot retrieve parameter from NULL component");

double result(0.0);
if(const auto group = boost::dynamic_pointer_cast<const Geometry::DetectorGroup>(comp))
{
const auto dets = group->getDetectors();
double avg(0.0);
for(auto it = dets.begin(); it!= dets.end(); ++it)
{
auto param = pmap.getRecursive((*it)->getComponentID(), name);
if(param) avg += param->value<double>();
else throw std::invalid_argument("ComptonProfile - Unable to find DetectorGroup component parameter \"" + name + "\".");
}
result = avg/static_cast<double>(group->nDets());
}
else
{
auto param = pmap.getRecursive(comp->getComponentID(), name);
if(param)
{
result = param->value<double>();
}
else
{
throw std::invalid_argument("ComptonProfile - Unable to find component parameter \"" + name + "\".");
}
}
return result;
}

/**
*/
ComptonProfile::ComptonProfile() : API::ParamFunction(), API::IFunction1D(),
Expand Down Expand Up @@ -89,18 +130,19 @@ namespace CurveFitting
}

DetectorParams detpar;
const auto & pmap = workspace->constInstrumentParameters();
detpar.l1 = sample->getDistance(*source);
detpar.l2 = det->getDistance(*sample);
detpar.theta = workspace->detectorTwoTheta(det);
detpar.t0 = getComponentParameter(*det,"t0")*1e-6; // Convert to seconds
detpar.efixed = getComponentParameter(*det,"efixed");
detpar.t0 = getComponentParameter(det, pmap, "t0")*1e-6; // Convert to seconds
detpar.efixed = getComponentParameter(det, pmap, "efixed");

ResolutionParams respar;
respar.dl1 = getComponentParameter(*det, "sigma_l1");
respar.dl2 = getComponentParameter(*det, "sigma_l2");
respar.dthe = getComponentParameter(*det,"sigma_theta"); //radians
respar.dEnLorentz = getComponentParameter(*det, "hwhm_lorentz");
respar.dEnGauss = getComponentParameter(*det, "sigma_gauss");
respar.dl1 = getComponentParameter(det, pmap, "sigma_l1");
respar.dl2 = getComponentParameter(det, pmap, "sigma_l2");
respar.dthe = getComponentParameter(det, pmap, "sigma_theta"); //radians
respar.dEnLorentz = getComponentParameter(det, pmap, "hwhm_lorentz");
respar.dEnGauss = getComponentParameter(det, pmap, "sigma_gauss");

this->cacheYSpaceValues(workspace->readX(m_wsIndex), workspace->isHistogramData(), detpar, respar);
}
Expand Down Expand Up @@ -304,25 +346,5 @@ namespace CurveFitting
std::transform(voigt.begin(), voigt.end(), voigt.begin(), std::bind2nd(std::multiplies<double>(), norm));
}

/**
* @param comp A reference to the component that should contain the parameter
* @param name The name of the parameter
* @returns The value of the parameter if it exists
* @throws A std::invalid_argument error if the parameter does not exist
*/
double ComptonProfile::getComponentParameter(const Geometry::IComponent & comp,const std::string &name) const
{
std::vector<double> pars = comp.getNumberParameter(name);
if(!pars.empty())
{
return pars[0];
}
else
{
throw std::invalid_argument("ComptonProfile - Unable to find component parameter \"" + name + "\".");
}
}


} // namespace CurveFitting
} // namespace Mantid

0 comments on commit 1745c65

Please sign in to comment.