Refs #6970 better number of point for profile

Code/Mantid/Framework/MDAlgorithms/src/IntegratePeaksMD.cpp

@@ -72,6 +72,9 @@ IntegratePeaksMD(InputWorkspace='TOPAZ_3131_md', PeaksWorkspace='peaks',
 #include "MantidMDAlgorithms/IntegratePeaksMD.h"
 #include "MantidMDEvents/CoordTransformDistance.h"
 #include "MantidKernel/ListValidator.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidAPI/AnalysisDataService.h"
 
 namespace Mantid
 {
@@ -193,13 +196,25 @@ namespace MDAlgorithms
     double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
     /// Cylinder Length to use around peaks for cylinder
     double cylinderLength = getProperty("CylinderLength");
+    Workspace2D_sptr ws2D;
+    size_t numSteps;
+    double deltaQ;
+    bool cylinderBool = getProperty("Cylinder");
+    if (cylinderBool)
+    {
+        numSteps = 15;
+        deltaQ = cylinderLength/static_cast<double>(numSteps-1);
+        size_t histogramNumber = peakWS->getNumberPeaks();
+        Workspace_sptr wsFit= WorkspaceFactory::Instance().create("Workspace2D",histogramNumber,numSteps,numSteps);
+        ws2D = boost::dynamic_pointer_cast<Workspace2D>(wsFit);
+        AnalysisDataService::Instance().addOrReplace("ProfilesToFit", ws2D);
+    }
     double backgroundCylinder = cylinderLength;
     double percentBackground = getProperty("PercentBackground");
     cylinderLength *= 1.0 - (percentBackground/100.);
     /// Replace intensity with 0
     bool replaceIntensity = getProperty("ReplaceIntensity");
     bool integrateEdge = getProperty("IntegrateIfOnEdge");
-    bool cylinderBool = getProperty("Cylinder");
     if (BackgroundInnerRadius < PeakRadius)
       BackgroundInnerRadius = PeakRadius;
 
@@ -319,28 +334,25 @@ namespace MDAlgorithms
 			// Perform the integration into whatever box is contained within.
 			std::vector<signal_t> signal_fit;
 
-			int numSteps = 20;
-			double deltaQ = cylinderLength/static_cast<double>(numSteps-1);
 			signal_fit.clear();
-			for (int j=0; j<numSteps; j++)signal_fit.push_back(0.0);
+			for (size_t j=0; j<numSteps; j++)signal_fit.push_back(0.0);
 			ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(PeakRadius), static_cast<coord_t>(cylinderLength), signal, errorSquared, signal_fit);
-			for (size_t j=0; j<signal_fit.size(); j++)std::cout << signal_fit[j]<<"  ";
-			std::cout <<i<<"\n";
-
 
 			// Integrate around the background radius
 			if (BackgroundOuterRadius > PeakRadius || percentBackground > 0.0)
 			{
 				// Get the total signal inside "BackgroundOuterRadius"
 
 				if (BackgroundOuterRadius < PeakRadius ) BackgroundOuterRadius = PeakRadius;
-				int numSteps=static_cast<int>((backgroundCylinder/deltaQ) + 1);
 				signal_fit.clear();
-				for (int j=0; j<numSteps; j++)signal_fit.push_back(0.0);
+				for (size_t j=0; j<numSteps; j++)signal_fit.push_back(0.0);
 				ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundOuterRadius), static_cast<coord_t>(backgroundCylinder), bgSignal, bgErrorSquared, signal_fit);
-				for (size_t j=0; j<signal_fit.size(); j++)std::cout << signal_fit[j]<<"  ";
-				std::cout <<i<<"\n";
-
+				for (size_t j = 0; j < numSteps; j++)
+				{
+					 ws2D->dataX(i)[j] = -0.5*backgroundCylinder + static_cast<double>(j) * deltaQ;
+					 ws2D->dataY(i)[j] = signal_fit[j];
+					 ws2D->dataE(i)[j] = std::sqrt(signal_fit[j]);
+				}
 
 				// Evaluate the signal inside "BackgroundInnerRadius"
 				signal_t interiorSignal = 0;
@@ -350,8 +362,6 @@ namespace MDAlgorithms
 				if (BackgroundInnerRadius != PeakRadius)
 				{
 					ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundInnerRadius), static_cast<coord_t>(cylinderLength), interiorSignal, interiorErrorSquared, signal_fit);
-					for (size_t j=0; j<signal_fit.size(); j++)std::cout << signal_fit[j]<<"  ";
-					std::cout <<i<<"\n";
 				}
 				else
 				{
@@ -381,6 +391,15 @@ namespace MDAlgorithms
 		        // But we add the errors together
 		        errorSquared += bgErrorSquared;
 			}
+			else
+			{
+				for (size_t j = 0; j < numSteps; j++)
+				{
+					 ws2D->dataX(i)[j] = -0.5*cylinderLength + static_cast<double>(j) * deltaQ;
+					 ws2D->dataY(i)[j] = signal_fit[j];
+					 ws2D->dataE(i)[j] = std::sqrt(signal_fit[j]);
+				}
+			}
       }
 
 

Code/Mantid/Framework/MDEvents/src/MDBox.cpp

@@ -622,7 +622,7 @@ namespace MDEvents
     const std::vector<MDE> & events = this->getConstEvents();
     typename std::vector<MDE>::const_iterator it = events.begin();
     typename std::vector<MDE>::const_iterator it_end = events.end();
-    int numSteps = 100;
+    size_t numSteps = signal_fit.size();
     double deltaQ = length/static_cast<double>(numSteps-1);
 
     // For each MDLeanEvent
@@ -633,7 +633,7 @@ namespace MDEvents
       if (out[0] < radius && std::fabs(out[1]) < 0.5*length)
       {
         // add event to appropriate y channel
-        int xchannel = static_cast<int>(floor((out[1] / deltaQ))+0.5) + (numSteps / 2);
+        size_t xchannel = static_cast<int>(floor((out[1] / deltaQ))+0.5) + (numSteps / 2);
         if (xchannel >= 0 || xchannel < numSteps ) signal_fit[xchannel] += static_cast<signal_t>(it->getSignal());
 
         signal += static_cast<signal_t>(it->getSignal());

Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp

@@ -1354,7 +1354,7 @@ namespace MDEvents
     size_t boxIndex[nd]; Utils::NestedForLoop::SetUp(nd, boxIndex, 0);
     size_t indexMaker[nd]; Utils::NestedForLoop::SetUpIndexMaker(nd, indexMaker, split);
 
-    int numSteps = 100;
+    size_t numSteps = signal_fit.size();
     double deltaQ = length/static_cast<double>(numSteps-1);
     bool allDone = false;
     while (!allDone)
@@ -1423,7 +1423,7 @@ namespace MDEvents
         coord_t out[nd];
         radiusTransform.apply(boxCenter, out);
         // add event to appropriate y channel
-        int xchannel = static_cast<int>(floor((out[1] / deltaQ))+0.5) + (numSteps / 2);
+        size_t xchannel = static_cast<int>(floor((out[1] / deltaQ))+0.5) + (numSteps / 2);
         if (xchannel >= 0 || xchannel < numSteps ) signal_fit[xchannel] += box->getSignal();
         // Use the integrated sum of signal in the box
         signal += box->getSignal();