Refs #9445. PoldiCalculateSpectrum2D now integrates peak profiles

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

@@ -7,6 +7,7 @@
 #include "MantidAPI/MultiDomainFunction.h"
 #include "MantidDataObjects/TableWorkspace.h"
 #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
+#include "MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h"
 
 
 namespace Mantid
@@ -49,13 +50,21 @@ class MANTID_SINQ_DLL PoldiCalculateSpectrum2D  : public API::Algorithm
     virtual int version() const;
     virtual const std::string category() const;
 
+protected:
+    PoldiPeakCollection_sptr getPeakCollection(TableWorkspace_sptr peakTable);
+    PoldiPeakCollection_sptr getIntegratedPeakCollection(PoldiPeakCollection_sptr rawPeakCollection);
+    PoldiPeakCollection_sptr getNormalizedPeakCollection(PoldiPeakCollection_sptr peakCollection);
+
+    boost::shared_ptr<MultiDomainFunction> getMultiDomainFunctionFromPeakCollection(PoldiPeakCollection_sptr peakCollection);
+
+    PoldiTimeTransformer_sptr m_timeTransformer;
+    double m_deltaT;
+
 private:
     virtual void initDocs();
     void init();
     void exec();
 
-    PoldiPeakCollection_sptr getPeakCollection(TableWorkspace_sptr peakTable);
-    boost::shared_ptr<MultiDomainFunction> getMultiDomainFunctionFromPeakCollection(PoldiPeakCollection_sptr peakCollection);
 };
 
 

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

@@ -11,39 +11,45 @@ TODO: Enter a full wiki-markup description of your algorithm here. You can then
 #include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
 
 #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
+#include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
+#include "MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h"
+#include "MantidAPI/IPeakFunction.h"
 
 namespace Mantid
 {
 namespace Poldi
 {
   // Register the algorithm into the AlgorithmFactory
   DECLARE_ALGORITHM(PoldiCalculateSpectrum2D)
-  
+
   using namespace API;
 
 
   //----------------------------------------------------------------------------------------------
   /** Constructor
    */
-  PoldiCalculateSpectrum2D::PoldiCalculateSpectrum2D()
+  PoldiCalculateSpectrum2D::PoldiCalculateSpectrum2D():
+      Algorithm(),
+      m_timeTransformer(),
+      m_deltaT(0.0)
   {
   }
-    
+
   //----------------------------------------------------------------------------------------------
   /** Destructor
    */
   PoldiCalculateSpectrum2D::~PoldiCalculateSpectrum2D()
   {
   }
-  
+
 
   //----------------------------------------------------------------------------------------------
   /// Algorithm's name for identification. @see Algorithm::name
   const std::string PoldiCalculateSpectrum2D::name() const { return "PoldiCalculateSpectrum2D";}
-  
+
   /// Algorithm's version for identification. @see Algorithm::version
   int PoldiCalculateSpectrum2D::version() const { return 1;}
-  
+
   /// Algorithm's category for identification. @see Algorithm::category
   const std::string PoldiCalculateSpectrum2D::category() const { return "SINQ\\Poldi\\PoldiSet";}
 
@@ -89,20 +95,33 @@ namespace Poldi
   void PoldiCalculateSpectrum2D::exec()
   {
       TableWorkspace_sptr peakTable = getProperty("PoldiPeakWorkspace");
-
       if(!peakTable) {
           throw std::runtime_error("Cannot proceed without peak workspace.");
       }
 
+      MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
+      m_deltaT = ws->readX(0)[1] - ws->readX(0)[0];
+
+      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));
+
       PoldiPeakCollection_sptr peakCollection = getPeakCollection(peakTable);
 
       if(peakCollection->intensityType() != PoldiPeakCollection::Integral) {
-          throw std::runtime_error("Peaks intensities need to be integral.");
+          /* Intensities must be integral for fitting procedure, so they have to be integrated
+           * and put into a new PoldiPeakCollection.
+           */
+
+          peakCollection = getIntegratedPeakCollection(peakCollection);
       }
 
-      boost::shared_ptr<MultiDomainFunction> mdFunction = getMultiDomainFunctionFromPeakCollection(peakCollection);
+      peakCollection = getNormalizedPeakCollection(peakCollection);
 
-      MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
+      boost::shared_ptr<MultiDomainFunction> mdFunction = getMultiDomainFunctionFromPeakCollection(peakCollection);
 
       size_t nspec = ws->getNumberHistograms();
       std::vector<size_t> wsi(nspec);
@@ -137,8 +156,74 @@ namespace Poldi
           throw std::runtime_error("Could not initialize peak collection.");
       }
   }
-  // TODO Auto-generated execute stub
 
+  PoldiPeakCollection_sptr PoldiCalculateSpectrum2D::getIntegratedPeakCollection(PoldiPeakCollection_sptr rawPeakCollection)
+  {
+      /* If no profile function is specified, it's not possible to get integrated
+       * intensities at all and we need to abort at this point.
+       */
+      if(!rawPeakCollection->hasProfileFunctionName()) {
+          throw std::runtime_error("Cannot integrate peak profiles without profile function.");
+      }
+
+      PeakFunctionIntegrator peakIntegrator(1e-10);
+
+      PoldiPeakCollection_sptr integratedPeakCollection(new PoldiPeakCollection(PoldiPeakCollection::Integral));
+      integratedPeakCollection->setProfileFunctionName(rawPeakCollection->getProfileFunctionName());
+
+      for(size_t i = 0; i < rawPeakCollection->peakCount(); ++i) {
+          PoldiPeak_sptr peak = rawPeakCollection->peak(i);
+
+          /* The integration is performed in time dimension,
+           * so the fwhm needs to be transformed.
+           */
+          double fwhmTime = m_timeTransformer->dToTOF(peak->fwhm(PoldiPeak::AbsoluteD));
+
+          IPeakFunction_sptr profileFunction = boost::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction(rawPeakCollection->getProfileFunctionName()));
+          profileFunction->setHeight(peak->intensity());
+          profileFunction->setFwhm(fwhmTime);
+
+          /* Because the integration is running from -inf to inf, it is necessary
+           * to set the centre to 0. Otherwise the transformation performed by
+           * the integration routine will create problems.
+           */
+          profileFunction->setCentre(0.0);
+
+          IntegrationResult integration = peakIntegrator.integrateInfinity(profileFunction);
+
+          if(!integration.success) {
+              throw std::runtime_error("Problem during peak integration. Aborting.");
+          }
+
+          PoldiPeak_sptr integratedPeak = peak->clone();
+          /* Integration is performed in the time domain and later everything is normalized
+           * by deltaT. In the original code this is done at this point, so this behavior is kept
+           * for now.
+           */
+          integratedPeak->setIntensity(UncertainValue(integration.result / m_deltaT));
+          integratedPeakCollection->addPeak(integratedPeak);
+      }
+
+      return integratedPeakCollection;
+  }
+
+  PoldiPeakCollection_sptr PoldiCalculateSpectrum2D::getNormalizedPeakCollection(PoldiPeakCollection_sptr peakCollection)
+  {
+      PoldiPeakCollection_sptr normalizedPeakCollection(new PoldiPeakCollection(PoldiPeakCollection::Integral));
+      normalizedPeakCollection->setProfileFunctionName(peakCollection->getProfileFunctionName());
+
+      for(size_t i = 0; i < peakCollection->peakCount(); ++i) {
+          PoldiPeak_sptr peak = peakCollection->peak(i);
+          double calculatedIntensity = m_timeTransformer->calculatedTotalIntensity(peak->d());
+
+          PoldiPeak_sptr normalizedPeak = peak->clone();
+          normalizedPeak->setIntensity(peak->intensity() / calculatedIntensity);
+
+          normalizedPeakCollection->addPeak(normalizedPeak);
+      }
+
+      return normalizedPeakCollection;
+  }
 
 
 } // namespace Poldi