Refs #9445. Refactored PoldiCalculateSpectrum2D.

Most parts of the algorithm are tested by the unit tests now. The actual spectrum calculation will be done in a system test,
as instrument definitions etc. are required for this.

Code/Mantid/Framework/SINQ/CMakeLists.txt

@@ -93,6 +93,7 @@ set ( TEST_FILES
 	Poldi2DFunctionTest.h
 	PoldiAutoCorrelationCoreTest.h
 	PoldiBasicChopperTest.h
+        PoldiCalculateSpectrum2DTest.h
 	PoldiChopperFactoryTest.h
         PoldiConversionsTest.h
 	PoldiDeadWireDecoratorTest.h

Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiCalculateSpectrum2D.h

@@ -58,6 +58,18 @@ class MANTID_SINQ_DLL PoldiCalculateSpectrum2D  : public API::Algorithm
 
     boost::shared_ptr<Poldi2DFunction> getFunctionFromPeakCollection(PoldiPeakCollection_sptr peakCollection);
 
+    API::MatrixWorkspace_sptr calculateSpectrum(PoldiPeakCollection_sptr peakCollection, API::MatrixWorkspace_sptr matrixWorkspace);
+
+    void setTimeTransformerFromInstrument(PoldiInstrumentAdapter_sptr poldiInstrument);
+    void setTimeTransformer(PoldiTimeTransformer_sptr poldiTimeTransformer);
+
+    void setDeltaTFromWorkspace(API::MatrixWorkspace_sptr matrixWorkspace);
+    void setDeltaT(double newDeltaT);
+    bool isValidDeltaT(double deltaT);
+
+    void throwOnInsufficientState();
+
+
     PoldiTimeTransformer_sptr m_timeTransformer;
     double m_deltaT;
 

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

@@ -8,13 +8,18 @@
 #include "MantidSINQ/PoldiUtilities/PoldiAbstractChopper.h"
 #include "MantidSINQ/PoldiUtilities/PoldiSourceSpectrum.h"
 #include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
+#include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
 
 #include "MantidSINQ/PoldiUtilities/PoldiHeliumDetector.h"
 #include "MantidSINQ/PoldiUtilities/PoldiConversions.h"
 #include "MantidGeometry/Instrument.h"
 #include "MantidGeometry/Instrument/FitParameter.h"
 #include "MantidKernel/Interpolation.h"
 
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidDataObjects/TableWorkspace.h"
+#include "MantidAPI/TableRow.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
 
 using ::testing::Return;
 
@@ -84,10 +89,10 @@ class  MockChopper : public PoldiAbstractChopper
 public:
     MockChopper() : PoldiAbstractChopper()
     {
-        double slits [] = {0.000000, 0.162156};
+        double slits [] = {0.000000, 0.162156, 0.250867, 0.3704, 0.439811, 0.588455, 0.761389, 0.895667};
         m_slitPositions = std::vector<double>(slits, slits + sizeof(slits) / sizeof(slits[0]));
 
-        double times [] = {0.000000, 243.234};
+        double times [] = {0.000000, 243.234, 376.3, 555.6, 659.716, 882.682, 1142.08, 1343.5};
         m_slitTimes = std::vector<double>(times, times + sizeof(times) / sizeof(times[0]));
     }
 
@@ -408,6 +413,112 @@ class FakePoldiInstrumentAdapter : public PoldiInstrumentAdapter
     }
 };
 
+
+class PoldiPeakCollectionHelpers
+{
+/* This class contains some static helper function to create
+ * peak collections and related components for testing purposes.
+ */
+public:
+
+/**
+ * This function creates a TableWorkspace which can be used by PoldiPeakCollection
+ *
+ * A TableWorkspace with four peaks is generated. The data comes from a run of the
+ * original analysis software on poldi2013n006904 (available in system tests and
+ * usage tests). Only the first four peaks are considered and no errors are provided.
+ *
+ * @return TableWorkspace in the format PoldiPeakCollection requires
+ */
+static DataObjects::TableWorkspace_sptr createPoldiPeakTableWorkspace()
+{
+    DataObjects::TableWorkspace_sptr tableWs = boost::dynamic_pointer_cast<DataObjects::TableWorkspace>(API::WorkspaceFactory::Instance().createTable());
+    tableWs->addColumn("str", "HKL");
+    tableWs->addColumn("str", "d");
+    tableWs->addColumn("str", "Q");
+    tableWs->addColumn("str", "Intensity");
+    tableWs->addColumn("str", "FWHM (rel.)");
+
+    tableWs->logs()->addProperty<std::string>("IntensityType", "Maximum");
+    tableWs->logs()->addProperty<std::string>("ProfileFunctionName", "Gaussian");
+
+    API::TableRow newRow = tableWs->appendRow();
+    newRow << "0 0 0" << "1.108644" << "5.667449" << "3286.152" << "0.002475747";
+
+    newRow = tableWs->appendRow();
+    newRow << "0 0 0" << "1.637539" << "3.836968" << "2951.696" << "0.002516417";
+
+    newRow = tableWs->appendRow();
+    newRow << "0 0 0" << "1.920200" << "3.272152" << "3238.473" << "0.002444439";
+
+    newRow = tableWs->appendRow();
+    newRow << "0 0 0" << "1.245958" << "5.042856" << "2219.592" << "0.002696334";
+
+    return tableWs;
+}
+
+/**
+ * This function creates a PoldiPeakCollection
+ *
+ * A PoldiPeakCollection is created with the same information as in createPoldiPeakTableWorkspace().
+ *
+ * @return PoldiPeakCollection with four example peaks.
+ */
+static PoldiPeakCollection_sptr createPoldiPeakCollectionMaximum()
+{
+    PoldiPeakCollection_sptr peaks(new PoldiPeakCollection);
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.108644), UncertainValue(3286.152), UncertainValue(0.002475747)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.637539), UncertainValue(2951.696), UncertainValue(0.002516417)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.920200), UncertainValue(3238.473), UncertainValue(0.002444439)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.245958), UncertainValue(2219.592), UncertainValue(0.002696334)));
+
+    peaks->setProfileFunctionName("Gaussian");
+
+    return peaks;
+}
+
+/**
+ * This function creates a PoldiPeakCollection with integrated intensities
+ *
+ * The same peaks as above, with integrated intensities in "channel units".
+ *
+ * @return PoldiPeakCollection with four example peaks.
+ */
+static PoldiPeakCollection_sptr createPoldiPeakCollectionIntegral()
+{
+    PoldiPeakCollection_sptr peaks(new PoldiPeakCollection(PoldiPeakCollection::Integral));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.108644), UncertainValue(15835.28906), UncertainValue(0.002475747)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.637539), UncertainValue(21354.32226), UncertainValue(0.002516417)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.920200), UncertainValue(26687.36132), UncertainValue(0.002444439)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.245958), UncertainValue(13091.51855), UncertainValue(0.002696334)));
+
+    peaks->setProfileFunctionName("Gaussian");
+
+    return peaks;
+}
+
+/**
+ * This function creates a PoldiPeakCollection with normalized intensities
+ *
+ * Again, the same peaks as above, but with normalized intensities (for 8 chopper slits)
+ *
+ * @return PoldiPeakCollection with four example peaks.
+ */
+static PoldiPeakCollection_sptr createPoldiPeakCollectionNormalized()
+{
+    PoldiPeakCollection_sptr peaks(new PoldiPeakCollection(PoldiPeakCollection::Integral));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.108644), UncertainValue(1.926395655), UncertainValue(0.002475747)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.920200), UncertainValue(4.773980141), UncertainValue(0.002444439)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.637539), UncertainValue(9.370919228), UncertainValue(0.002516417)));
+    peaks->addPeak(PoldiPeak::create(MillerIndices(0, 0, 0), UncertainValue(1.245958), UncertainValue(1.758037806), UncertainValue(0.002696334)));
+
+    peaks->setProfileFunctionName("Gaussian");
+
+    return peaks;
+}
+
+};
+
 }
 }
 #endif // POLDIMOCKINSTRUMENTHELPERS_H

Code/Mantid/Framework/SINQ/src/PoldiCalculateSpectrum2D.cpp

@@ -6,6 +6,7 @@ TODO: Enter a full wiki-markup description of your algorithm here. You can then
 
 #include "MantidDataObjects/Workspace2D.h"
 #include "MantidDataObjects/TableWorkspace.h"
+#include "MantidAPI/TableRow.h"
 #include "MantidAPI/FunctionFactory.h"
 #include "MantidAPI/MultiDomainFunction.h"
 #include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
@@ -93,6 +94,11 @@ namespace Poldi
       return mdFunction;
   }
 
+  void PoldiCalculateSpectrum2D::setTimeTransformer(PoldiTimeTransformer_sptr poldiTimeTransformer)
+  {
+      m_timeTransformer = poldiTimeTransformer;
+  }
+
   void PoldiCalculateSpectrum2D::exec()
   {
       TableWorkspace_sptr peakTable = getProperty("PoldiPeakWorkspace");
@@ -101,36 +107,73 @@ namespace Poldi
       }
 
       MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
-      m_deltaT = ws->readX(0)[1] - ws->readX(0)[0];
+      setDeltaTFromWorkspace(ws);
 
-      if(m_deltaT <= 0.0) {
-          throw std::runtime_error("Invalid time bin size.");
-      }
-
-      PoldiInstrumentAdapter_sptr poldiInstrument(new PoldiInstrumentAdapter(ws->getInstrument(), ws->run()));
-      m_timeTransformer = PoldiTimeTransformer_sptr(new PoldiTimeTransformer(poldiInstrument));
+      setTimeTransformerFromInstrument(PoldiInstrumentAdapter_sptr(new PoldiInstrumentAdapter(ws)));
 
       PoldiPeakCollection_sptr peakCollection = getPeakCollection(peakTable);
 
-      if(peakCollection->intensityType() != PoldiPeakCollection::Integral) {
-          /* Intensities must be integral for fitting procedure, so they have to be integrated
-           * and put into a new PoldiPeakCollection.
-           */
+      setProperty("OutputWorkspace", calculateSpectrum(peakCollection, ws));
+  }
 
-          peakCollection = getIntegratedPeakCollection(peakCollection);
-      }
+  MatrixWorkspace_sptr PoldiCalculateSpectrum2D::calculateSpectrum(PoldiPeakCollection_sptr peakCollection, MatrixWorkspace_sptr matrixWorkspace)
+  {
+      PoldiPeakCollection_sptr integratedPeaks = getIntegratedPeakCollection(peakCollection);
+      PoldiPeakCollection_sptr normalizedPeakCollection = getNormalizedPeakCollection(integratedPeaks);
 
-      peakCollection = getNormalizedPeakCollection(peakCollection);
+      boost::shared_ptr<IFunction> mdFunction = getFunctionFromPeakCollection(normalizedPeakCollection);
+      IAlgorithm_sptr fit = createChildAlgorithm("Fit", -1, -1, true);
 
-      boost::shared_ptr<IFunction> mdFunction = getFunctionFromPeakCollection(peakCollection);
+      if(!fit) {
+          throw std::runtime_error("Could not initialize 'Fit'-algorithm.");
+      }
 
-      IAlgorithm_sptr fit = createChildAlgorithm("Fit", -1, -1, true);
       fit->setProperty("Function", boost::dynamic_pointer_cast<IFunction>(mdFunction));
-      fit->setProperty("InputWorkspace", ws);
+      fit->setProperty("InputWorkspace", matrixWorkspace);
       fit->setProperty("CreateOutput", true);
       fit->setProperty("MaxIterations", 0);
+      fit->setProperty("Minimizer", "Levenberg-MarquardtMD");
 
       fit->execute();
+
+      MatrixWorkspace_sptr outputWs = fit->getProperty("OutputWorkspace");
+
+      return outputWs;
+  }
+
+  void PoldiCalculateSpectrum2D::setTimeTransformerFromInstrument(PoldiInstrumentAdapter_sptr poldiInstrument)
+  {
+      setTimeTransformer(PoldiTimeTransformer_sptr(new PoldiTimeTransformer(poldiInstrument)));
+  }
+
+  void PoldiCalculateSpectrum2D::setDeltaTFromWorkspace(MatrixWorkspace_sptr matrixWorkspace)
+  {
+      if(matrixWorkspace->getNumberHistograms() < 1) {
+          throw std::invalid_argument("MatrixWorkspace does not contain any data.");
+      }
+
+      MantidVec xData = matrixWorkspace->readX(0);
+
+      if(xData.size() < 2) {
+          throw std::invalid_argument("Cannot process MatrixWorkspace with less than 2 x-values.");
+      }
+
+      // difference between first and second x-value is assumed to be the bin width.
+      setDeltaT(matrixWorkspace->readX(0)[1] - matrixWorkspace->readX(0)[0]);
+  }
+
+  void PoldiCalculateSpectrum2D::setDeltaT(double newDeltaT)
+  {
+      if(!isValidDeltaT(newDeltaT)) {
+          throw std::invalid_argument("Time bin size must be larger than 0.");
+      }
+
+      m_deltaT = newDeltaT;
+  }
+
+  bool PoldiCalculateSpectrum2D::isValidDeltaT(double deltaT)
+  {
+      return deltaT > 0.0;
   }
 
   PoldiPeakCollection_sptr PoldiCalculateSpectrum2D::getPeakCollection(TableWorkspace_sptr peakTable)
@@ -144,6 +187,26 @@ namespace Poldi
 
   PoldiPeakCollection_sptr PoldiCalculateSpectrum2D::getIntegratedPeakCollection(PoldiPeakCollection_sptr rawPeakCollection)
   {
+      if(!rawPeakCollection) {
+          throw std::invalid_argument("Cannot proceed with invalid PoldiPeakCollection.");
+      }
+
+      if(!isValidDeltaT(m_deltaT)) {
+          throw std::invalid_argument("Cannot proceed with invalid time bin size.");
+      }
+
+      if(!m_timeTransformer) {
+          throw std::invalid_argument("Cannot proceed with invalid PoldiTimeTransformer.");
+      }
+
+      if(rawPeakCollection->intensityType() == PoldiPeakCollection::Integral) {
+          /* Intensities are integral already - don't need to do anything,
+           * except cloning the collection, to make behavior consistent, since
+           * integrating also results in a new peak collection.
+           */
+          return rawPeakCollection->clone();
+      }
+
       /* If no profile function is specified, it's not possible to get integrated
        * intensities at all and we need to abort at this point.
        */
@@ -185,6 +248,7 @@ namespace Poldi
            * by deltaT. In the original code this is done at this point, so this behavior is kept
            * for now.
            */
+          integratedPeak->setD(integratedPeak->d() - 0.0005);
           integratedPeak->setIntensity(UncertainValue(integration.result / m_deltaT));
           integratedPeakCollection->addPeak(integratedPeak);
       }
@@ -194,6 +258,14 @@ namespace Poldi
 
   PoldiPeakCollection_sptr PoldiCalculateSpectrum2D::getNormalizedPeakCollection(PoldiPeakCollection_sptr peakCollection)
   {
+      if(!peakCollection) {
+          throw std::invalid_argument("Cannot proceed with invalid PoldiPeakCollection.");
+      }
+
+      if(!m_timeTransformer) {
+          throw std::invalid_argument("Cannot proceed without PoldiTimeTransformer.");
+      }
+
       PoldiPeakCollection_sptr normalizedPeakCollection(new PoldiPeakCollection(PoldiPeakCollection::Integral));
       normalizedPeakCollection->setProfileFunctionName(peakCollection->getProfileFunctionName());