From 06606cc87c57f610512cb29b61537f57ef9ebdd8 Mon Sep 17 00:00:00 2001 From: Martyn Gigg Date: Wed, 11 Dec 2013 15:35:36 +0000 Subject: [PATCH] Refactoring to avoid code duplication. * getComponentParameter is now only defined in ConvertToYSpace * DetectorParams struct has moved to ConvertToYSpace * Y calculation has moved to ConvertToYSpace Refs #8598 --- .../CalculateGammaBackground.h | 7 +- .../inc/MantidCurveFitting/ComptonProfile.h | 41 +----- .../ComptonScatteringCountRate.h | 4 +- .../inc/MantidCurveFitting/ConvertToYSpace.h | 23 +++- .../src/CalculateGammaBackground.cpp | 38 +++--- .../CurveFitting/src/ComptonProfile.cpp | 119 +++++----------- .../src/ComptonScatteringCountRate.cpp | 2 + .../CurveFitting/src/ConvertToYSpace.cpp | 127 +++++++++++++++++- .../src/GramCharlierComptonProfile.cpp | 2 + .../CurveFitting/test/ConvertToYSpaceTest.h | 13 ++ .../WorkspaceCreationHelper.h | 2 +- .../src/WorkspaceCreationHelper.cpp | 7 +- 12 files changed, 227 insertions(+), 158 deletions(-) diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h index 9bd2a19e11af..52e0affe6af1 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateGammaBackground.h @@ -3,14 +3,17 @@ #include "MantidAPI/Algorithm.h" -#include "MantidCurveFitting/ComptonProfile.h" - #include namespace Mantid { namespace CurveFitting { + //--------------------------------------------------------------------------- + // Forward declarations + //--------------------------------------------------------------------------- + struct DetectorParams; + struct ResolutionParams; /** diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonProfile.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonProfile.h index fc076b667f8a..9c4c8acb9b1f 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonProfile.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonProfile.h @@ -10,17 +10,12 @@ namespace Mantid { namespace CurveFitting { - /// Simple data structure to store nominal detector values - /// It avoids some functions taking a huge number of arguments - struct DetectorParams - { - double l1; ///< source-sample distance in metres - double l2; ///< sample-detector distance in metres - double theta; ///< scattering angle in radians - double t0; ///< time delay in microseconds - double efixed; ///< final energy - }; + //--------------------------------------------------------------------------- + // Forward declarations + //--------------------------------------------------------------------------- + struct DetectorParams; + //--------------------------------------------------------------------------- /// Simple data structure to store resolution parameter values /// It avoids some functions taking a huge number of arguments struct ResolutionParams @@ -58,10 +53,6 @@ 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(); @@ -125,28 +116,6 @@ namespace CurveFitting /// Store the mass values double m_mass; - /// Source to sample distance - double m_l1; - /// Std dev of l1 distance - double m_sigmaL1; - /// Sample to detector distance - double m_l2; - /// Std dev of l1 distance - double m_sigmaL2; - /// Theta value for this spectrum - double m_theta; - /// Std Dev theta value for this spectrum - double m_sigmaTheta; - - /// Final energy - double m_e1; - /// T0 value for this spectrum in seconds - double m_t0; - /// Lorentzian HWHM of the foil analyser energy - double m_hwhmGaussE; - /// Gaussian HWHM of the foil analyser energy - double m_hwhmLorentzE; - /// Voigt function boost::shared_ptr m_voigt; diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonScatteringCountRate.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonScatteringCountRate.h index 8e2df7b7b4dc..34bfda722528 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonScatteringCountRate.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonScatteringCountRate.h @@ -67,9 +67,9 @@ namespace CurveFitting void cacheBackground(const API::IFunction1D_sptr & profile, const size_t paramsOffset); /// Set up the constraint matrices - void createConstraintMatrices(const MantidVec & xValues); + void createConstraintMatrices(const std::vector & xValues); /// Set up positivity constraint matrix - void createPositivityCM(const MantidVec & xValues); + void createPositivityCM(const std::vector & xValues); /// Set up equality constraint matrix void createEqualityCM(const size_t nmasses); diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h index 3a115034d7d7..3fe3949578f5 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h @@ -7,10 +7,17 @@ namespace Mantid { namespace CurveFitting { - //--------------------------------------------------------------------------- - // Forward declarations - //--------------------------------------------------------------------------- - struct DetectorParams; + + /// Simple data structure to store nominal detector values + /// It avoids some functions taking a huge number of arguments + struct DetectorParams + { + double l1; ///< source-sample distance in metres + double l2; ///< sample-detector distance in metres + double theta; ///< scattering angle in radians + double t0; ///< time delay in microseconds + double efixed; ///< final energy + }; /** Takes a workspace with X axis in TOF and converts it to Y-space where the transformation is defined @@ -46,6 +53,10 @@ namespace CurveFitting int version() const; const std::string category() const; + /// Retrieve a component parameter + static double getComponentParameter(const Geometry::IComponent_const_sptr & comp, + const Geometry::ParameterMap &pmap, + const std::string &name); /// Convert single time value to Y,Q & Ei values static void calculateY(double & yspace, double & qspace, double &ei, const double mass, const double tmicro, @@ -57,6 +68,8 @@ namespace CurveFitting void init(); void exec(); + /// Perform the conversion to Y-space + void convert(const size_t i); /// Check and store appropriate input data void retrieveInputs(); /// Create the output workspace @@ -66,6 +79,8 @@ namespace CurveFitting /// Input workspace API::MatrixWorkspace_sptr m_inputWS; + /// The input mass in AMU + double m_mass; /// Source-sample distance double m_l1; /// Sample position diff --git a/Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp b/Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp index abb94f850f09..dd884c7f7428 100644 --- a/Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/CalculateGammaBackground.cpp @@ -4,6 +4,8 @@ *WIKI*/ #include "MantidCurveFitting/CalculateGammaBackground.h" +#include "MantidCurveFitting/ComptonProfile.h" +#include "MantidCurveFitting/ConvertToYSpace.h" #include "MantidAPI/CompositeFunction.h" #include "MantidAPI/FunctionProperty.h" @@ -226,15 +228,15 @@ namespace Mantid detPar.l1 = m_l1; detPar.l2 = m_samplePos.distance(detPos); detPar.theta = m_inputWS->detectorTwoTheta(det); //radians - detPar.t0 = ComptonProfile::getComponentParameter(det, pmap, "t0")*1e-6; // seconds - detPar.efixed = ComptonProfile::getComponentParameter(det, pmap,"efixed"); + detPar.t0 = ConvertToYSpace::getComponentParameter(det, pmap, "t0")*1e-6; // seconds + detPar.efixed = ConvertToYSpace::getComponentParameter(det, pmap,"efixed"); ResolutionParams detRes; - 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"); + detRes.dl1 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l1"); // DL0 + detRes.dl2 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l2"); // DL1 + detRes.dthe = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_theta"); //DTH in radians + detRes.dEnGauss = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_gauss"); + detRes.dEnLorentz = ConvertToYSpace::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 @@ -264,15 +266,15 @@ namespace Mantid detPar.l1 = m_l1; detPar.l2 = m_samplePos.distance(detPos); detPar.theta = m_inputWS->detectorTwoTheta(det); //radians - detPar.t0 = ComptonProfile::getComponentParameter(det, pmap, "t0")*1e-6; // seconds - detPar.efixed = ComptonProfile::getComponentParameter(det, pmap,"efixed"); + detPar.t0 = ConvertToYSpace::getComponentParameter(det, pmap, "t0")*1e-6; // seconds + detPar.efixed = ConvertToYSpace::getComponentParameter(det, pmap,"efixed"); ResolutionParams detRes; - 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"); + detRes.dl1 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l1"); // DL0 + detRes.dl2 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l2"); // DL1 + detRes.dthe = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_theta"); //DTH in radians + detRes.dEnGauss = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_gauss"); + detRes.dEnLorentz = ConvertToYSpace::getComponentParameter(det, pmap, "hwhm_lorentz"); const size_t nxvalues = m_backgroundWS->blocksize(); std::vector foilSpectrum(nxvalues); @@ -537,16 +539,16 @@ namespace Mantid FoilInfo descr; descr.thetaMin = thetaRng0.first; descr.thetaMax = thetaRng0.second; - descr.lorentzWidth = ComptonProfile::getComponentParameter(foil0, pmap, "hwhm_lorentz"); - descr.gaussWidth = ComptonProfile::getComponentParameter(foil0, pmap, "sigma_gauss"); + descr.lorentzWidth = ConvertToYSpace::getComponentParameter(foil0, pmap, "hwhm_lorentz"); + descr.gaussWidth = ConvertToYSpace::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 = ComptonProfile::getComponentParameter(foil1, pmap, "hwhm_lorentz"); - descr.gaussWidth = ComptonProfile::getComponentParameter(foil1, pmap, "sigma_gauss"); + descr.lorentzWidth = ConvertToYSpace::getComponentParameter(foil1, pmap, "hwhm_lorentz"); + descr.gaussWidth = ConvertToYSpace::getComponentParameter(foil1, pmap, "sigma_gauss"); m_foils1[i] = descr; //copy } diff --git a/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp b/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp index 698e7fefbd32..6b63700ee92e 100644 --- a/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp @@ -2,8 +2,8 @@ // Includes //----------------------------------------------------------------------------- #include "MantidCurveFitting/ComptonProfile.h" +#include "MantidCurveFitting/ConvertToYSpace.h" #include "MantidAPI/FunctionFactory.h" -#include "MantidGeometry/Instrument/DetectorGroup.h" #include namespace Mantid @@ -21,54 +21,11 @@ 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(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(); - else throw std::invalid_argument("ComptonProfile - Unable to find DetectorGroup component parameter \"" + name + "\"."); - } - result = avg/static_cast(group->nDets()); - } - else - { - auto param = pmap.getRecursive(comp->getComponentID(), name); - if(param) - { - result = param->value(); - } - else - { - throw std::invalid_argument("ComptonProfile - Unable to find component parameter \"" + name + "\"."); - } - } - return result; - } - /** */ ComptonProfile::ComptonProfile() : API::ParamFunction(), API::IFunction1D(), m_log(Kernel::Logger::get("ComptonProfile")), - m_wsIndex(0), m_mass(0.0), m_l1(0.0), m_sigmaL1(0.0), - m_l2(0.0), m_sigmaL2(0.0), m_theta(0.0), m_sigmaTheta(0.0), m_e1(0.0), - m_t0(0.0), m_hwhmGaussE(0.0), m_hwhmLorentzE(0.0), m_voigt(), + m_wsIndex(0), m_mass(0.0), m_voigt(), m_yspace(), m_modQ(), m_e0(),m_resolutionSigma(0.0), m_lorentzFWHM(0.0) {} @@ -134,15 +91,15 @@ namespace CurveFitting detpar.l1 = sample->getDistance(*source); detpar.l2 = det->getDistance(*sample); detpar.theta = workspace->detectorTwoTheta(det); - detpar.t0 = getComponentParameter(det, pmap, "t0")*1e-6; // Convert to seconds - detpar.efixed = getComponentParameter(det, pmap, "efixed"); + detpar.t0 = ConvertToYSpace::getComponentParameter(det, pmap, "t0")*1e-6; // Convert to seconds + detpar.efixed = ConvertToYSpace::getComponentParameter(det, pmap, "efixed"); ResolutionParams respar; - 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"); + respar.dl1 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l1"); + respar.dl2 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l2"); + respar.dthe = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_theta"); //radians + respar.dEnLorentz = ConvertToYSpace::getComponentParameter(det, pmap, "hwhm_lorentz"); + respar.dEnGauss = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_gauss"); this->cacheYSpaceValues(workspace->readX(m_wsIndex), workspace->isHistogramData(), detpar, respar); } @@ -156,27 +113,19 @@ namespace CurveFitting void ComptonProfile::cacheYSpaceValues(const std::vector & tseconds, const bool isHistogram, const DetectorParams & detpar, const ResolutionParams & respar) { - //geometry - m_l1 = detpar.l1; - m_l2 = detpar.l2; - m_theta = detpar.theta; - m_e1 = detpar.efixed; - m_t0 = detpar.t0; - //resolution - m_sigmaL1 = respar.dl1; - m_sigmaL2 = respar.dl2; - m_sigmaTheta = respar.dthe; - m_hwhmLorentzE = respar.dEnLorentz; - m_hwhmGaussE = STDDEV_TO_HWHM*respar.dEnGauss; + // geometry + double theta = detpar.theta; //cache for frequent access + double hwhmLorentzE = respar.dEnLorentz; + double hwhmGaussE = STDDEV_TO_HWHM*respar.dEnGauss; // ------ Fixed coefficients related to resolution & Y-space transforms ------------------ const double mn = PhysicalConstants::NeutronMassAMU; const double mevToK = PhysicalConstants::E_mev_toNeutronWavenumberSq; const double massToMeV = 0.5*PhysicalConstants::NeutronMass/PhysicalConstants::meV; // Includes factor of 1/2 - const double v1 = std::sqrt(m_e1/massToMeV); - const double k1 = std::sqrt(m_e1/mevToK); - const double l2l1 = m_l2/m_l1; + const double v1 = std::sqrt(detpar.efixed/massToMeV); + const double k1 = std::sqrt(detpar.efixed/mevToK); + const double l2l1 = detpar.l2/detpar.l1; // Resolution dependence @@ -186,13 +135,13 @@ namespace CurveFitting if((m_mass-1.0) > DBL_EPSILON) { double x0(0.0),x1(0.0); - gsl_poly_solve_quadratic(m_mass-1.0, 2.0*std::cos(m_theta), -(m_mass+1.0), &x0, &x1); + gsl_poly_solve_quadratic(m_mass-1.0, 2.0*std::cos(theta), -(m_mass+1.0), &x0, &x1); k0k1 = std::max(x0,x1); // K0/K1 at y=0 } else { // solution is simply s = 1/cos(theta) - k0k1 = 1.0/std::cos(m_theta); + k0k1 = 1.0/std::cos(theta); } double qy0(0.0), wgauss(0.0); @@ -201,25 +150,25 @@ namespace CurveFitting qy0 = std::sqrt(k1*k1*m_mass*(k0k1*k0k1 - 1)); double k0k1p3 = std::pow(k0k1,3); double r1 = -(1.0 + l2l1*k0k1p3); - double r2 = 1.0 - l2l1*k0k1p3 + l2l1*std::pow(k0k1,2)*std::cos(m_theta) - k0k1*std::cos(m_theta); + double r2 = 1.0 - l2l1*k0k1p3 + l2l1*std::pow(k0k1,2)*std::cos(theta) - k0k1*std::cos(theta); double factor = (0.2413/qy0)*((m_mass/mn)*r1 - r2); - m_lorentzFWHM = std::abs(factor*m_hwhmLorentzE*2); - wgauss = std::abs(factor*m_hwhmGaussE*2); + m_lorentzFWHM = std::abs(factor*hwhmLorentzE*2); + wgauss = std::abs(factor*hwhmGaussE*2); } else { - qy0 = k1*std::tan(m_theta); - double factor = (0.2413*2.0/k1)*std::abs((std::cos(m_theta) + l2l1)/std::sin(m_theta)); - m_lorentzFWHM = m_hwhmLorentzE*factor; - wgauss = m_hwhmGaussE*factor; + qy0 = k1*std::tan(theta); + double factor = (0.2413*2.0/k1)*std::abs((std::cos(theta) + l2l1)/std::sin(theta)); + m_lorentzFWHM = hwhmLorentzE*factor; + wgauss = hwhmGaussE*factor; } double k0y0 = k1*k0k1; // k0_y0 = k0 value at y=0 - double wtheta = 2.0*STDDEV_TO_HWHM*std::abs(k0y0*k1*std::sin(m_theta)/qy0)*m_sigmaTheta; - double common = (m_mass/mn) - 1 + k1*std::cos(m_theta)/k0y0; - double wl1 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,2)/(qy0*m_l1))*common)*m_sigmaL1; - double wl2 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,3)/(k1*qy0*m_l1))*common)*m_sigmaL2; + double wtheta = 2.0*STDDEV_TO_HWHM*std::abs(k0y0*k1*std::sin(theta)/qy0)*respar.dthe; + double common = (m_mass/mn) - 1 + k1*std::cos(theta)/k0y0; + double wl1 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,2)/(qy0*detpar.l1))*common)*respar.dl1; + double wl2 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,3)/(k1*qy0*detpar.l1))*common)*respar.dl2; m_resolutionSigma = std::sqrt(std::pow(wgauss,2) + std::pow(wtheta,2) + std::pow(wl1,2) + std::pow(wl2,2)); @@ -239,14 +188,8 @@ namespace CurveFitting for(size_t i = 0; i < nData; ++i) { const double tsec = (isHistogram) ? 0.5*(tseconds[i] + tseconds[i+1]) : tseconds[i]; - const double v0 = m_l1/(tsec - m_t0 - (m_l2/v1)); - const double ei = massToMeV*v0*v0; - m_e0[i] = ei; - const double w = ei - m_e1; - const double k0 = std::sqrt(ei/mevToK); - const double q = std::sqrt(k0*k0 + k1*k1 - 2.0*k0*k1*std::cos(m_theta)); - m_modQ[i] = q; - m_yspace[i] = 0.2393*(m_mass/q)*(w - (mevToK*q*q/m_mass)); + ConvertToYSpace::calculateY(m_yspace[i], m_modQ[i],m_e0[i], + m_mass,tsec*1e06,k1,v1,detpar); } } diff --git a/Code/Mantid/Framework/CurveFitting/src/ComptonScatteringCountRate.cpp b/Code/Mantid/Framework/CurveFitting/src/ComptonScatteringCountRate.cpp index a36f2b5ac822..e2b2c5272755 100644 --- a/Code/Mantid/Framework/CurveFitting/src/ComptonScatteringCountRate.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/ComptonScatteringCountRate.cpp @@ -5,6 +5,8 @@ #include +#include + namespace Mantid { namespace CurveFitting diff --git a/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp b/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp index 30dce83c50ed..16bf6bcaca94 100644 --- a/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp @@ -3,9 +3,9 @@ TODO: Enter a full wiki-markup description of your algorithm here. You can then *WIKI*/ #include "MantidCurveFitting/ConvertToYSpace.h" -#include "MantidCurveFitting/ComptonProfile.h" #include "MantidAPI/WorkspaceValidators.h" +#include "MantidGeometry/Instrument/DetectorGroup.h" #include "MantidKernel/BoundedValidator.h" #include @@ -31,7 +31,7 @@ namespace CurveFitting /** Constructor */ ConvertToYSpace::ConvertToYSpace() - : Algorithm(), m_inputWS(), m_l1(0.0), m_samplePos(),m_outputWS() + : Algorithm(), m_inputWS(), m_mass(0.0), m_l1(0.0), m_samplePos(),m_outputWS() { } @@ -50,11 +50,57 @@ namespace CurveFitting /// Sets documentation strings for this algorithm void ConvertToYSpace::initDocs() { - this->setWikiSummary("Converts in workspace in units of TOF to Y-space as defined in Compton scattering field"); - this->setOptionalMessage("Converts in workspace in units of TOF to Y-space as defined in Compton scattering field"); + this->setWikiSummary("Converts workspace in units of TOF to Y-space as defined in Compton scattering field"); + this->setOptionalMessage("Converts workspace in units of TOF to Y-space as defined in Compton scattering field"); } //---------------------------------------------------------------------------------------------- + + /** + * 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 ConvertToYSpace::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(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(); + else + throw std::invalid_argument("ComptonProfile - Unable to find DetectorGroup component parameter \"" + name + "\"."); + } + result = avg/static_cast(group->nDets()); + } + else + { + auto param = pmap.getRecursive(comp->getComponentID(), name); + if(param) + { + result = param->value(); + } + else + { + throw std::invalid_argument("ComptonProfile - Unable to find component parameter \"" + name + "\"."); + } + } + return result; + } + + //---------------------------------------------------------------------------------------------- + /** * @param yspace Output yspace value * @param qspace Output qspace value @@ -89,8 +135,9 @@ namespace CurveFitting void ConvertToYSpace::init() { auto wsValidator = boost::make_shared(); - wsValidator->add("TOF"); + wsValidator->add(false); // point data wsValidator->add(); + wsValidator->add("TOF"); declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input,wsValidator), "An input workspace."); @@ -102,6 +149,9 @@ namespace CurveFitting declareProperty(new WorkspaceProperty<>("OutputWorkspace","",Direction::Output), "An output workspace."); } + + + //---------------------------------------------------------------------------------------------- /** Execute the algorithm. */ @@ -109,14 +159,81 @@ namespace CurveFitting { retrieveInputs(); createOutputWorkspace(); + + const int64_t nhist = static_cast(m_inputWS->getNumberHistograms()); + const int64_t nreports = nhist; + auto progress = boost::make_shared(this, 0.0, 1.0, nreports); + + PARALLEL_FOR2(m_inputWS, m_outputWS) + for(int64_t i = 0; i < nhist; ++i) + { + PARALLEL_START_INTERUPT_REGION + + try + { + convert(i); + } + catch(Exception::NotFoundError &) + { + g_log.warning("No detector defined for index=" + boost::lexical_cast(i) + ". Zeroing spectrum."); + m_outputWS->maskWorkspaceIndex(i); + } + + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + + setProperty("OutputWorkspace", m_outputWS); + } + + /** + * Convert the spectrum at the given index on the input workspace + * and place the output in the pre-allocated output workspace + * @param index Index on the input & output workspaces giving the spectrum to convert + */ + void ConvertToYSpace::convert(const size_t index) + { + auto det = m_inputWS->getDetector(index); + 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, pmap, "t0"); // micro-seconds + detPar.efixed = getComponentParameter(det, pmap,"efixed"); + const double v1 = std::sqrt(detPar.efixed/MASS_TO_MEV); + const double k1 = std::sqrt(detPar.efixed/PhysicalConstants::E_mev_toNeutronWavenumberSq); + + auto & outX = m_outputWS->dataX(index); + auto & outY = m_outputWS->dataY(index); + auto & outE = m_outputWS->dataE(index); + const auto & inX = m_inputWS->dataX(index); + const auto & inY = m_inputWS->dataY(index); + const auto & inE = m_inputWS->dataE(index); + + const size_t npts = inY.size(); + for(size_t j = 0; j < npts; ++j) + { + double ys(0.0),qs(0.0),ei(0.0); + calculateY(ys,qs,ei,m_mass,inX[j],k1,v1,detPar); + outX[j] = ys; + const double prefactor = qs/pow(ei,0.1); + outY[j] = prefactor*inY[j]; + outE[j] = prefactor*inE[j]; + } } + /** * Caches input details for the peak information */ void ConvertToYSpace::retrieveInputs() { m_inputWS = getProperty("InputWorkspace"); + m_mass = getProperty("Mass"); cacheInstrumentGeometry(); } diff --git a/Code/Mantid/Framework/CurveFitting/src/GramCharlierComptonProfile.cpp b/Code/Mantid/Framework/CurveFitting/src/GramCharlierComptonProfile.cpp index 4d10384c24a4..24c8c287d473 100644 --- a/Code/Mantid/Framework/CurveFitting/src/GramCharlierComptonProfile.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/GramCharlierComptonProfile.cpp @@ -9,6 +9,8 @@ #include // for factorial #include +#include + namespace Mantid { namespace CurveFitting diff --git a/Code/Mantid/Framework/CurveFitting/test/ConvertToYSpaceTest.h b/Code/Mantid/Framework/CurveFitting/test/ConvertToYSpaceTest.h index ea8f702b0a64..80d8153741fb 100644 --- a/Code/Mantid/Framework/CurveFitting/test/ConvertToYSpaceTest.h +++ b/Code/Mantid/Framework/CurveFitting/test/ConvertToYSpaceTest.h @@ -62,6 +62,19 @@ class ConvertToYSpaceTest : public CxxTest::TestSuite TS_ASSERT_THROWS(alg->setProperty("InputWorkspace", testWS), std::invalid_argument); } + void test_Input_Workspace_In_TOF_With_Instrument_But_No_Detector_Parameters_Throws_Error_On_Execution() + { + auto alg = createAlgorithm(); + auto testWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(1,10); + testWS->getAxis(0)->setUnit("TOF"); + + alg->setProperty("InputWorkspace", testWS); + alg->setProperty("Mass", 1.0097); + alg->setRethrows(true); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } + private: Mantid::API::IAlgorithm_sptr createAlgorithm() { diff --git a/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h b/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h index d59aee9ed81d..e0ecbbe8e4b1 100644 --- a/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h +++ b/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h @@ -162,7 +162,7 @@ namespace WorkspaceCreationHelper * Data filled with: Y: 2.0, E: sqrt(2.0), X: nbins of width 1 starting at 0 */ Mantid::DataObjects::Workspace2D_sptr create2DWorkspaceWithFullInstrument(int nHist, int nBins, - bool includeMonitors = false, bool startYNegative = false); + bool includeMonitors = false, bool startYNegative = false, bool isHistogram = false); /** * Create a test workspace with a Theta numeric axis instead of a spectrum axis diff --git a/Code/Mantid/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp b/Code/Mantid/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp index 1adea60151ed..5529e4fcaeb8 100644 --- a/Code/Mantid/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp +++ b/Code/Mantid/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp @@ -312,14 +312,17 @@ namespace WorkspaceCreationHelper * pervious. * Data filled with: Y: 2.0, E: sqrt(2.0), X: nbins of width 1 starting at 0 */ - Workspace2D_sptr create2DWorkspaceWithFullInstrument(int nhist, int nbins, bool includeMonitors, bool startYNegative) + Workspace2D_sptr create2DWorkspaceWithFullInstrument(int nhist, int nbins, bool includeMonitors, + bool startYNegative, bool isHistogram) { if( includeMonitors && nhist < 2 ) { throw std::invalid_argument("Attemping to 2 include monitors for a workspace with fewer than 2 histograms"); } - Workspace2D_sptr space = Create2DWorkspaceBinned(nhist, nbins); // A 1:1 spectra is created by default + Workspace2D_sptr space; + if(isHistogram) space = Create2DWorkspaceBinned(nhist, nbins); // A 1:1 spectra is created by default + else space = Create2DWorkspace123(nhist,nbins, false); space->setTitle("Test histogram"); // actually adds a property call run_title to the logs space->getAxis(0)->setUnit("TOF"); space->setYUnit("Counts");