Refs #9445. Added PoldiSpectrumDomainFunction class

Code/Mantid/Framework/SINQ/CMakeLists.txt

@@ -24,6 +24,7 @@ set ( SRC_FILES
         src/PoldiUtilities/PoldiInstrumentAdapter.cpp
 	src/PoldiUtilities/PoldiPeak.cpp
 	src/PoldiUtilities/PoldiPeakCollection.cpp
+        src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp
 	src/PoldiUtilities/UncertainValue.cpp
 	src/ProjectMD.cpp
 	src/SINQHMListener.cpp
@@ -62,6 +63,7 @@ set ( INC_FILES
 	inc/MantidSINQ/PoldiUtilities/PoldiMockInstrumentHelpers.h
 	inc/MantidSINQ/PoldiUtilities/PoldiPeak.h
 	inc/MantidSINQ/PoldiUtilities/PoldiPeakCollection.h
+        inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h
 	inc/MantidSINQ/PoldiUtilities/UncertainValue.h
 	inc/MantidSINQ/PoldiUtilities/UncertainValueIO.h
 	inc/MantidSINQ/ProjectMD.h

Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h

@@ -0,0 +1,140 @@
+#ifndef MANTID_SINQ_POLDISPECTRUMDOMAINFUNCTION_H_
+#define MANTID_SINQ_POLDISPECTRUMDOMAINFUNCTION_H_
+
+#include "MantidSINQ/DllConfig.h"
+#include "MantidAPI/ParamFunction.h"
+#include <string>
+
+#include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
+#include "MantidSINQ/PoldiUtilities/PoldiSourceSpectrum.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+#include "MantidDataObjects/Workspace2D.h"
+
+namespace Mantid
+{
+namespace Poldi
+{
+
+/** PoldiSpectrumDomainFunction : TODO: DESCRIPTION
+
+    Copyright &copy; 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+    This file is part of Mantid.
+
+    Mantid is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.
+
+    Mantid is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+    File change history is stored at: <https://github.com/mantidproject/mantid>
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+  */
+using namespace API;
+using namespace DataObjects;
+
+struct DetectorElementCharacteristics {
+    DetectorElementCharacteristics()
+    { }
+
+    DetectorElementCharacteristics(int element, PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper)
+    {
+        distance = detector->distanceFromSample(element);
+        totalDistance = detector->distanceFromSample(element) + chopper->distanceFromSample();
+        twoTheta = detector->twoTheta(element);
+        sinTheta = sin(twoTheta / 2.0);
+        cosTheta = cos(twoTheta / 2.0);
+        tof1A = 4947.990234375;//2.0 * totalDistance * sinTheta * PhysicalConstants::NeutronMass / (PhysicalConstants::h * 1e7);//4947.99013618171464192258; //063732507753820628;//4947.990; //Conversions::dtoTOF(1.0, totalDistance, sinTheta);
+    }
+
+    double distance;
+    double totalDistance;
+    double twoTheta;
+    double sinTheta;
+    double cosTheta;
+    double tof1A;
+};
+
+class DetectorElementData
+{
+public:
+    DetectorElementData(int element, DetectorElementCharacteristics center, PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper)
+    {
+        DetectorElementCharacteristics current(element, detector, chopper);
+
+        m_intensityFactor = pow(center.distance / current.distance, 2.0) * current.sinTheta / center.sinTheta;
+        m_lambdaFactor = 2.0 * current.sinTheta / center.tof1A;
+        m_timeFactor = current.sinTheta / center.sinTheta * current.totalDistance / center.totalDistance;
+        m_widthFactor = current.cosTheta - center.cosTheta;
+    }
+
+    double intensityFactor() const { return m_intensityFactor; }
+    double lambdaFactor() const { return m_lambdaFactor; }
+    double timeFactor() const { return m_timeFactor; }
+    double widthFactor() const { return m_widthFactor; }
+
+protected:
+    double m_intensityFactor;
+    double m_lambdaFactor;
+    double m_timeFactor;
+    double m_widthFactor;
+};
+
+typedef boost::shared_ptr<const DetectorElementData> DetectorElementData_const_sptr;
+
+class MANTID_SINQ_DLL PoldiSpectrumDomainFunction : public ParamFunction
+{
+public:
+    PoldiSpectrumDomainFunction();
+    virtual ~PoldiSpectrumDomainFunction()
+    {}
+
+    virtual std::string name() const { return "PoldiSpectrumDomainFunction"; }
+
+    virtual void setWorkspace(boost::shared_ptr<const Workspace> ws);
+    virtual void function(const FunctionDomain &domain, FunctionValues &values) const;
+
+
+protected:
+    virtual void init();
+
+    double getArrivalTime(double tof);
+
+    void initializeParametersFromWorkspace(Workspace2D_const_sptr workspace2D);
+    void initializeFromInstrument(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper);
+
+    std::vector<double> getChopperSlitOffsets(PoldiAbstractChopper_sptr chopper);
+    std::vector<DetectorElementData_const_sptr> getDetectorElementData(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper);
+    DetectorElementCharacteristics getDetectorCenterCharacteristics(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper);
+
+    double dToTOF(double d) const;
+    double timeTransformedWidth(double widthT, size_t detectorIndex) const;
+    double timeTransformedCentre(double centreT, size_t detectorIndex) const;
+    double timeTransformedIntensity(double areaD, double centreT, size_t detectorIndex) const;
+    double detectorElementIntensity(double centreT, size_t detectorIndex) const;
+
+    double actualFunction(double x, double x0, double sigma, double area) const;
+
+    PoldiSourceSpectrum_const_sptr m_spectrum;
+    DetectorElementCharacteristics m_detectorCenter;
+    std::vector<DetectorElementData_const_sptr> m_detectorElementData;
+
+    std::vector<double> m_chopperSlitOffsets;
+
+    double m_detectorEfficiency;
+    double m_deltaT;
+};
+
+
+} // namespace Poldi
+} // namespace Mantid
+
+#endif  /* MANTID_SINQ_POLDISPECTRUMDOMAINFUNCTION_H_ */

Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp

@@ -0,0 +1,190 @@
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
+
+#include "MantidAPI/FunctionFactory.h"
+#include "MantidAPI/Workspace.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include <stdexcept>
+
+#include "MantidAPI/FunctionDomain1D.h"
+
+namespace Mantid
+{
+namespace Poldi
+{
+
+using namespace DataObjects;
+using namespace API;
+
+DECLARE_FUNCTION(PoldiSpectrumDomainFunction)
+
+PoldiSpectrumDomainFunction::PoldiSpectrumDomainFunction() :
+    ParamFunction()
+{
+
+}
+
+void PoldiSpectrumDomainFunction::initializeParametersFromWorkspace(Workspace2D_const_sptr workspace2D)
+{
+    PoldiInstrumentAdapter adapter(workspace2D->getInstrument(),  workspace2D->run());
+    initializeFromInstrument(adapter.detector(), adapter.chopper());
+
+    m_spectrum = boost::const_pointer_cast<const PoldiSourceSpectrum>(adapter.spectrum());
+    m_deltaT = workspace2D->readX(0)[1] - workspace2D->readX(0)[1];
+}
+
+void PoldiSpectrumDomainFunction::initializeFromInstrument(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper)
+{
+    m_chopperSlitOffsets = getChopperSlitOffsets(chopper);
+    //m_chopperCycleTime = chopper->cycleTime();
+
+    m_detectorCenter = getDetectorCenterCharacteristics(detector, chopper);
+    m_detectorElementData = getDetectorElementData(detector, chopper);
+    m_detectorEfficiency = 0.88;
+}
+
+void PoldiSpectrumDomainFunction::setWorkspace(boost::shared_ptr<const Workspace> ws)
+{
+    Workspace2D_const_sptr workspace2D = boost::dynamic_pointer_cast<const Workspace2D>(ws);
+
+    if(!workspace2D) {
+        throw std::invalid_argument("PoldiSpectrumDomainFunction can only work with Workspace2D.");
+    }
+
+    initializeParametersFromWorkspace(workspace2D);
+}
+
+void PoldiSpectrumDomainFunction::function(const FunctionDomain &domain, FunctionValues &values) const
+{    
+    const FunctionDomain1DSpectrum &spectrumDomain = dynamic_cast<const FunctionDomain1DSpectrum &>(domain);
+
+    try{
+        spectrumDomain.getWorkspaceIndex();
+    } catch(...) {
+        throw std::invalid_argument("PoldiSpectrumDomainFunction can only work with FunctionDomain1DSpectrum.");
+    }
+
+    size_t index = spectrumDomain.getWorkspaceIndex();
+    int domainSize = static_cast<int>(spectrumDomain.size());
+
+    /* Parameters are given in d, but need to be processed in arrival time.
+     * This section performs the conversions. They depend on several factors
+     * terminated by the position in the detector, so the index is stored.
+     */
+    double fwhm = getParameter("Fwhm");
+    double fwhmT = timeTransformedWidth(dToTOF(fwhm), index);
+    double fwhmChannel = fwhmT / m_deltaT;
+    double sigmaChannel = fwhmChannel / (2.0 * sqrt(2.0 * log(2.0)));
+
+    double centre = getParameter("Centre");
+    double centreTRaw = dToTOF(centre);
+    double centreT = timeTransformedCentre(centreTRaw, index);
+
+    double area = getParameter("Area");
+    double areaT = timeTransformedIntensity(area, centreTRaw, index);
+
+    /* Once all the factors are all there, the calculation needs to be
+     * performed with one offset per chopper slit.
+     */
+    for(size_t o = 0; o < m_chopperSlitOffsets.size(); ++o) {
+        double centreTOffset = centreT + m_chopperSlitOffsets[o];
+        double centreTOffsetChannel = centreTOffset / m_deltaT;
+
+        /* Calculations are performed in channel units
+         * Needs to be signed integer, because the profile can extend beyond the left edge,
+         * which results in negative indices. Since the spectrum "wraps around" the
+         * indices have to be transformed to the right edge.
+         */
+        int centreChannel = static_cast<int>(centreTOffsetChannel);
+        int widthChannels = std::max(2, static_cast<int>(fwhmChannel * 2.0));
+
+        for(int i = centreChannel - widthChannels; i <= centreChannel + widthChannels; ++i) {
+            /* Since the POLDI spectra "wrap around" on the time axis, the x-value of the
+             * domain can not be used, because if the profile extends to x < 0, it should appear
+             * at 500 - x. The same for the other side.
+             */
+            int cleanChannel = i % domainSize;
+            if(cleanChannel < 0) {
+                cleanChannel += domainSize;
+            }
+
+            double xValue = static_cast<double>(i) + 0.5;
+
+            /* This is a workaround for later, when "actualFunction" will be replaced with
+             * an arbitrary profile.
+             */
+            values.addToCalculated(cleanChannel, actualFunction(xValue, centreTOffsetChannel, sigmaChannel, areaT));
+        }
+    }
+}
+
+void PoldiSpectrumDomainFunction::init() {
+    declareParameter("Area", 1.0);
+    declareParameter("Fwhm", 1.0);
+    declareParameter("Centre", 0.0);
+}
+
+std::vector<double> PoldiSpectrumDomainFunction::getChopperSlitOffsets(PoldiAbstractChopper_sptr chopper)
+{
+    const std::vector<double> &chopperSlitTimes = chopper->slitTimes();
+    std::vector<double> offsets;
+    offsets.reserve(chopperSlitTimes.size());
+    for(std::vector<double>::const_iterator time = chopperSlitTimes.begin(); time != chopperSlitTimes.end(); ++time) {
+        offsets.push_back(*time + chopper->zeroOffset());
+    }
+
+    return offsets;
+}
+
+std::vector<DetectorElementData_const_sptr> PoldiSpectrumDomainFunction::getDetectorElementData(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper)
+{
+    std::vector<DetectorElementData_const_sptr> data(detector->elementCount());
+
+    DetectorElementCharacteristics center = getDetectorCenterCharacteristics(detector, chopper);
+
+    for(int i = 0; i < static_cast<int>(detector->elementCount()); ++i) {
+        data[i] = DetectorElementData_const_sptr(new DetectorElementData(i, center, detector, chopper));
+    }
+
+    return data;
+}
+
+DetectorElementCharacteristics PoldiSpectrumDomainFunction::getDetectorCenterCharacteristics(PoldiAbstractDetector_sptr detector, PoldiAbstractChopper_sptr chopper)
+{
+    return DetectorElementCharacteristics(static_cast<int>(detector->centralElement()), detector, chopper);
+}
+
+double PoldiSpectrumDomainFunction::dToTOF(double d) const
+{
+    return m_detectorCenter.tof1A * d;
+}
+
+double PoldiSpectrumDomainFunction::timeTransformedWidth(double widthT, size_t detectorIndex) const
+{
+    return widthT;// + m_detectorElementData[detectorIndex]->widthFactor() * 0.0;
+}
+
+double PoldiSpectrumDomainFunction::timeTransformedCentre(double centreT, size_t detectorIndex) const
+{
+    return centreT * m_detectorElementData[detectorIndex]->timeFactor();
+}
+
+double PoldiSpectrumDomainFunction::timeTransformedIntensity(double areaD, double centreT, size_t detectorIndex) const
+{
+    return areaD * detectorElementIntensity(centreT, detectorIndex);
+}
+
+double PoldiSpectrumDomainFunction::detectorElementIntensity(double centreT, size_t detectorIndex) const
+{
+    double lambda = centreT * m_detectorElementData[detectorIndex]->lambdaFactor();
+    double intensity = m_spectrum->intensity(lambda) * m_detectorElementData[detectorIndex]->intensityFactor();
+
+    return intensity * (1.0 - exp(-m_detectorEfficiency * lambda));
+}
+
+double PoldiSpectrumDomainFunction::actualFunction(double x, double x0, double sigma, double area) const
+{
+    return area / (sqrt(2.0 * M_PI) * sigma) * exp(-0.5 * pow((x - x0) / sigma, 2.0));
+}
+
+} // namespace Poldi
+} // namespace Mantid