Refs #7449 rename FitPeakBackground and only output peak point range

Code/Mantid/Framework/Algorithms/CMakeLists.txt

@@ -103,6 +103,7 @@ set ( SRC_FILES
 	src/FindDeadDetectors.cpp
 	src/FindDetectorsOutsideLimits.cpp
 	src/FindPeaks.cpp
+	src/FindPeakBackground.cpp
 	src/FlatBackground.cpp
 	src/FlatPlateAbsorption.cpp
 	src/GeneralisedSecondDifference.cpp
@@ -184,7 +185,6 @@ set ( SRC_FILES
 	src/SassenaFFT.cpp
 	src/Scale.cpp
 	src/ScaleX.cpp
-	src/SeparateBackgroundFromSignal.cpp
 	src/SetUncertainties.cpp
 	src/ShiftLogTime.cpp
 	src/SignalOverError.cpp
@@ -325,6 +325,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/FindDeadDetectors.h
 	inc/MantidAlgorithms/FindDetectorsOutsideLimits.h
 	inc/MantidAlgorithms/FindPeaks.h
+	inc/MantidAlgorithms/FindPeakBackground.h
 	inc/MantidAlgorithms/FlatBackground.h
 	inc/MantidAlgorithms/FlatPlateAbsorption.h
 	inc/MantidAlgorithms/GSLFunctions.h
@@ -407,7 +408,6 @@ set ( INC_FILES
 	inc/MantidAlgorithms/SassenaFFT.h
 	inc/MantidAlgorithms/Scale.h
 	inc/MantidAlgorithms/ScaleX.h
-	inc/MantidAlgorithms/SeparateBackgroundFromSignal.h
 	inc/MantidAlgorithms/SetUncertainties.h
 	inc/MantidAlgorithms/ShiftLogTime.h
 	inc/MantidAlgorithms/SignalOverError.h
@@ -551,6 +551,7 @@ set ( TEST_FILES
 	FindDeadDetectorsTest.h
 	FindDetectorsOutsideLimitsTest.h
 	FindPeaksTest.h
+	FindPeakBackgroundTest.h
 	FlatBackgroundTest.h
 	FlatPlateAbsorptionTest.h
 	GenerateEventsFilterTest.h
@@ -622,7 +623,6 @@ set ( TEST_FILES
 	SassenaFFTTest.h
 	ScaleTest.h
 	ScaleXTest.h
-	SeparateBackgroundFromSignalTest.h
 	ShiftLogTimeTest.h
 	SignalOverErrorTest.h
 	SmoothDataTest.h

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SeparateBackgroundFromSignal.h → Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h

@@ -1,5 +1,5 @@
-#ifndef MANTID_ALGORITHMS_SeparateBackgroundFromSignal_H_
-#define MANTID_ALGORITHMS_SeparateBackgroundFromSignal_H_
+#ifndef MANTID_ALGORITHMS_FindPeakBackground_H_
+#define MANTID_ALGORITHMS_FindPeakBackground_H_
 
 #include "MantidKernel/System.h"
 #include "MantidAPI/Algorithm.h"
@@ -9,7 +9,7 @@ namespace Mantid
 namespace Algorithms
 {
 
-  /** SeparateBackgroundFromSignal : Calculate Zscore for a Matrix Workspace
+  /** FindPeakBackground : Calculate Zscore for a Matrix Workspace
     
     Copyright © 2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
 
@@ -31,14 +31,14 @@ namespace Algorithms
     File change history is stored at: <https://github.com/mantidproject/mantid>
     Code Documentation is available at: <http://doxygen.mantidproject.org>
   */
-  class DLLExport SeparateBackgroundFromSignal : public API::Algorithm
+  class DLLExport FindPeakBackground : public API::Algorithm
   {
   public:
-    SeparateBackgroundFromSignal();
-    virtual ~SeparateBackgroundFromSignal();
+    FindPeakBackground();
+    virtual ~FindPeakBackground();
 
     /// Algorithm's name for identification overriding a virtual method
-    virtual const std::string name() const { return "SeparateBackgroundFromSignal";}
+    virtual const std::string name() const { return "FindPeakBackground";}
 
     /// Algorithm's version for identification overriding a virtual method
     virtual int version() const { return 1;}
@@ -53,7 +53,7 @@ namespace Algorithms
     void init();
     /// Implement abstract Algorithm methods
     void exec();
-    double moment(MantidVec& X, size_t n, double mean, int k);
+    double moment4(MantidVec& X, size_t n, double mean);
 	struct cont_peak
 	{
 		size_t start;
@@ -74,4 +74,4 @@ namespace Algorithms
 } // namespace Algorithms
 } // namespace Mantid
 
-#endif  /* MANTID_ALGORITHMS_SeparateBackgroundFromSignal_H_ */
+#endif  /* MANTID_ALGORITHMS_FindPeakBackground_H_ */

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FindPeaks.h

@@ -91,6 +91,8 @@ class DLLExport FindPeaks : public API::Algorithm
   virtual int version() const { return (1); }
   /// Algorithm's category for identification
   virtual const std::string category() const { return "Optimization\\PeakFinding"; }
+  /// needed by FindPeaksBackground
+  int getVectorIndex(const MantidVec &vecX, double x);
 
 private:
   /// Sets documentation strings for this algorithm
@@ -122,9 +124,6 @@ class DLLExport FindPeaks : public API::Algorithm
   /// Fit peak
   void fitPeak(const API::MatrixWorkspace_sptr &input, const int spectrum, const int i_min, const int i_max, const int i_centre);
 
-
-  int getVectorIndex(const MantidVec &vecX, double x);
-
   void fitPeakHighBackground(const API::MatrixWorkspace_sptr &input, const size_t spectrum, const int &i_centre, const int &i_min, const int &i_max,
                              double &in_bg0, double &in_bg1, double &in_bg2);
 

Code/Mantid/Framework/Algorithms/src/SeparateBackgroundFromSignal.cpp → Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp

@@ -4,8 +4,8 @@ Separates background from signal for spectra of a workspace.
 
 *WIKI*/
 
-#include "MantidAlgorithms/SeparateBackgroundFromSignal.h"
-
+#include "MantidAlgorithms/FindPeakBackground.h"
+#include "MantidAlgorithms/FindPeaks.h"
 #include "MantidAPI/WorkspaceProperty.h"
 #include "MantidAPI/MatrixWorkspace.h"
 #include "MantidAPI/WorkspaceFactory.h"
@@ -26,26 +26,26 @@ namespace Mantid
 namespace Algorithms
 {
 
-  DECLARE_ALGORITHM(SeparateBackgroundFromSignal)
+  DECLARE_ALGORITHM(FindPeakBackground)
 
   //----------------------------------------------------------------------------------------------
   /** Constructor
    */
-  SeparateBackgroundFromSignal::SeparateBackgroundFromSignal()
+  FindPeakBackground::FindPeakBackground()
   {
   }
 
   //----------------------------------------------------------------------------------------------
   /** Destructor
    */
-  SeparateBackgroundFromSignal::~SeparateBackgroundFromSignal()
+  FindPeakBackground::~FindPeakBackground()
   {
   }
 
   //----------------------------------------------------------------------------------------------
   /** WIKI:
    */
-  void SeparateBackgroundFromSignal::initDocs()
+  void FindPeakBackground::initDocs()
   {
     setWikiSummary("Separates background from signal for spectra of a workspace.");
     setOptionalMessage("Separates background from signal for spectra of a workspace.");
@@ -54,27 +54,31 @@ namespace Algorithms
   //----------------------------------------------------------------------------------------------
   /** Define properties
     */
-  void SeparateBackgroundFromSignal::init()
+  void FindPeakBackground::init()
   {
     auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input);
-    declareProperty(inwsprop, "Name of input MatrixWorkspace to have Z-score calculated.");
-
-    auto outwsprop = new WorkspaceProperty<Workspace2D>("OutputWorkspace", "", Direction::Output);
-    declareProperty(outwsprop, "Name of the output Workspace2D containing the Z-scores.");
+    declareProperty(inwsprop, "Name of input MatrixWorkspace to have signal and background separated.");
 
-    declareProperty("WorkspaceIndex", EMPTY_INT(), "Index of the spectrum to have Z-score calculated. "
+    declareProperty("WorkspaceIndex", EMPTY_INT(), "Index of the spectrum to have peak and background separated. "
                     "Default is to calculate for all spectra.");
 
+    declareProperty(new ArrayProperty<double>("FitWindow"),
+                    "Optional: enter a comma-separated list of the expected X-position of window to fit. The number of values must be exactly two.");
+
+    auto outwsprop = new WorkspaceProperty<Workspace2D>("OutputWorkspace", "", Direction::Output);
+    declareProperty(outwsprop, "Name of the output Workspace2D containing the peak.");
+
   }
 
   //----------------------------------------------------------------------------------------------
   /** Execute body
     */
-  void SeparateBackgroundFromSignal::exec()
+  void FindPeakBackground::exec()
   {
     // 1. Get input and validate
     MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
     int inpwsindex = getProperty("WorkspaceIndex");
+    std::vector<double> m_vecFitWindows = getProperty("FitWindow");
 
     bool separateall = false;
     if (inpwsindex == EMPTY_INT())
@@ -94,13 +98,27 @@ namespace Algorithms
     }
     size_t sizex = inpWS->readX(0).size();
     size_t sizey = inpWS->readY(0).size();
+    size_t n = sizey;
+    size_t l0 = 0;
+    const MantidVec& inpX = inpWS->readX(0);
+
+    if (m_vecFitWindows.size() > 0)
+    {
+  	  Mantid::Algorithms::FindPeaks fp;
+  	  l0 = fp.getVectorIndex(inpX, m_vecFitWindows[0]);
+      n = fp.getVectorIndex(inpX, m_vecFitWindows[1]);
+      if (n < sizey) n++;
+    }
 
     Workspace2D_sptr outWS = boost::dynamic_pointer_cast<Workspace2D>(WorkspaceFactory::Instance().create(
-          "Workspace2D", numspec, sizex, sizey));
+          "Workspace2D", numspec, 2, 2));
 
-    // 3. Get Z values
+    // 3. Get Y values
+    Progress prog(this, 0, 1.0, numspec);
+    //PARALLEL_FOR2(inpWS, outWS)
     for (size_t i = 0; i < numspec; ++i)
     {
+      //PARALLEL_START_INTERUPT_REGION
       // a) figure out wsindex
       size_t wsindex;
       if (separateall)
@@ -128,94 +146,110 @@ namespace Algorithms
       MantidVec& vecX = outWS->dataX(i);
       MantidVec& vecY = outWS->dataY(i);
       MantidVec& vecE = outWS->dataE(i);
+
+	  // save output of mask * Y
+	  vecX[0] = inpX[l0];
+	  vecX[1] = inpX[n-1];
+	  vecY[0] = inpY[l0];
+	  vecY[1] = inpY[n-1];
+	  vecE[0] = inpE[l0];
+	  vecE[1] = inpE[n-1];
+
       double Ymean, Yvariance, Ysigma;
-      MantidVec maskedY = inpWS->readY(wsindex);
-      size_t n = maskedY.size();
-      MantidVec mask(n,0.0);
-      double xn = static_cast<double>(n);
+      MantidVec maskedY;
+      for (size_t l = l0; l < n; ++l)maskedY.push_back(inpY[l]);
+      MantidVec mask(n-l0,0.0);
+      double xn = static_cast<double>(n-l0);
       double k = 1.0;
       do
       {
     	  Statistics stats = getStatistics(maskedY);
     	  Ymean = stats.mean;
     	  Yvariance = stats.standard_deviation * stats.standard_deviation;
-    	  Ysigma = std::sqrt((moment(maskedY,n,Ymean,4)-(xn-3.0)/(xn-1.0) * moment(maskedY,n,Ymean,2))/xn);
+    	  Ysigma = std::sqrt((moment4(maskedY,n-l0,Ymean)-(xn-3.0)/(xn-1.0) * Yvariance)/xn);
 	      MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end());
 	      const size_t pos = it - maskedY.begin();
 		  maskedY[pos] = 0;
 		  mask[pos] = 1.0;
       }
       while (std::abs(Ymean-Yvariance) > k * Ysigma);
 
-      if(n > 5)
+      if(n-l0 > 5)
       {
     	  // remove single outliers
 		  if (mask[1] == mask[2] && mask[2] == mask[3])
 			  mask[0] = mask[1];
 		  if (mask[0] == mask[2] && mask[2] == mask[3])
 			  mask[1] = mask[2];
-		  for (size_t l = 2; l < n-3; ++l)
+		  for (size_t l = 2; l < n-l0-3; ++l)
 		  {
 			  if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2]))
 			  {
 				  mask[l] = mask[l+1];
 			  }
 		  }
-		  if (mask[n-2] == mask[n-3] && mask[n-3] == mask[n-4])
-			  mask[n-1] = mask[n-2];
-		  if (mask[n-1] == mask[n-3] && mask[n-3] == mask[n-4])
-			  mask[n-2] = mask[n-1];
+		  if (mask[n-l0-2] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
+			  mask[n-l0-1] = mask[n-l0-2];
+		  if (mask[n-l0-1] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
+			  mask[n-l0-2] = mask[n-l0-1];
 
 		  // mask regions not connected to largest region
 		  // for loop can start > 1 for multiple peaks
 		  vector<cont_peak> peaks;
-		  for (size_t l = 1; l < n; ++l)
+		  if (mask[0] == 1)
+		  {
+			  peaks.push_back(cont_peak());
+			  peaks[peaks.size()-1].start = l0;
+                  }
+		  for (size_t l = 1; l < n-l0; ++l)
 		  {
 			  if (mask[l] != mask[l-1] && mask[l] == 1)
 			  {
 				  peaks.push_back(cont_peak());
-				  peaks[peaks.size()-1].start = l;
+				  peaks[peaks.size()-1].start = l+l0;
 			  }
 			  if (peaks.size() > 0)
 			  {
 				  size_t ipeak = peaks.size()-1;
 				  if (mask[l] != mask[l-1] && mask[l] == 0)
 				  {
-					  peaks[ipeak].stop = l-1;
+					  peaks[ipeak].stop = l+l0-1;
 				  }
-				  if (inpY[l] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l];
+				  if (inpY[l+l0] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l+l0];
 			  }
 		  }
 		  if(peaks.size()> 0)
 		  {
 			  if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1;
 			  std::sort(peaks.begin(), peaks.end(), by_len());
-			  for (size_t l = 1; l < peaks.size(); ++l)
-			  {
-				  for (size_t j = peaks[l].start; j <= peaks[l].stop; ++j)
-				  {
-					  mask[j] = 0;
-				  }
-			  }
+
+			  // save endpoints
+			  vecX[0] = inpX[peaks[0].start];
+			  // extra point for histogram input
+			  vecX[1] = inpX[peaks[0].stop + sizex - sizey];
+			  vecY[0] = inpY[peaks[0].start];
+			  vecY[1] = inpY[peaks[0].stop];
+			  vecE[0] = inpE[peaks[0].start];
+			  vecE[1] = inpE[peaks[0].stop];
 		  }
       }
-	  // save output of mask * Y
-      vecX = inpX;
-      std::transform(mask.begin(), mask.end(), inpY.begin(), vecY.begin(), std::multiplies<double>());
-      std::transform(mask.begin(), mask.end(), inpE.begin(), vecE.begin(), std::multiplies<double>());
+
+	  prog.report();
+	  //PARALLEL_END_INTERUPT_REGION
     } // ENDFOR
+    //PARALLEL_CHECK_INTERUPT_REGION
 
     // 4. Set the output
     setProperty("OutputWorkspace", outWS);
 
     return;
   }
-  double SeparateBackgroundFromSignal::moment(MantidVec& X, size_t n, double mean, int k)
+  double FindPeakBackground::moment4(MantidVec& X, size_t n, double mean)
   {
 	  double sum=0.0;
 	  for (size_t i = 0; i < n; ++i)
 	  {
-		  sum += std::pow(X[i]-mean, k);
+		  sum += (X[i]-mean)*(X[i]-mean)*(X[i]-mean)*(X[i]-mean);
 	  }
 	  sum /= static_cast<double>(n);
           return sum;