Refs #6940 write profiles along cylinder

mantidproject · Jun 23, 2013 · bbd814e · bbd814e
1 parent a27e927
commit bbd814e
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 9 deletions.
diff --git a/Code/Mantid/Framework/MDAlgorithms/src/IntegratePeaksMD.cpp b/Code/Mantid/Framework/MDAlgorithms/src/IntegratePeaksMD.cpp
@@ -193,7 +193,6 @@ namespace MDAlgorithms
     double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
     /// Cylinder Length to use around peaks for cylinder
     double cylinderLength = getProperty("CylinderLength");
-    cylinderLength /= (2.0 * M_PI);
     double backgroundCylinder = cylinderLength;
     double percentBackground = getProperty("PercentBackground");
     cylinderLength *= 1.0 - (percentBackground/100.);
@@ -320,8 +319,8 @@ namespace MDAlgorithms
 			// Perform the integration into whatever box is contained within.
 			std::vector<signal_t> signal_fit;
 
-			double deltaQ = 0.004/(2*M_PI);
-			int numSteps=static_cast<int>((cylinderLength/deltaQ) + 1);
+			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);
 			ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(PeakRadius), static_cast<coord_t>(cylinderLength), signal, errorSquared, signal_fit);
@@ -335,7 +334,6 @@ namespace MDAlgorithms
 				// Get the total signal inside "BackgroundOuterRadius"
 
 				if (BackgroundOuterRadius < PeakRadius ) BackgroundOuterRadius = PeakRadius;
-				double deltaQ = 0.004/(2*M_PI);
 				int numSteps=static_cast<int>((backgroundCylinder/deltaQ) + 1);
 				signal_fit.clear();
 				for (int j=0; j<numSteps; j++)signal_fit.push_back(0.0);

diff --git a/Code/Mantid/Framework/MDEvents/src/MDBox.cpp b/Code/Mantid/Framework/MDEvents/src/MDBox.cpp
@@ -622,8 +622,8 @@ 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();
-    double deltaQ = 0.004/(2*M_PI);
-    int numSteps=static_cast<int>((length/deltaQ) + 1);
+    int numSteps = 100;
+    double deltaQ = length/static_cast<double>(numSteps-1);
 
     // For each MDLeanEvent
     for (; it != it_end; ++it)
@@ -632,10 +632,9 @@ namespace MDEvents
       radiusTransform.apply(it->getCenter(), out);
       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);
-        if (xchannel >= 0 || xchannel < numSteps ) signal_fit[xchannel]++;
+        if (xchannel >= 0 || xchannel < numSteps ) signal_fit[xchannel] += static_cast<signal_t>(it->getSignal());
 
         signal += static_cast<signal_t>(it->getSignal());
         errorSquared += static_cast<signal_t>(it->getErrorSquared());

diff --git a/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp b/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp
@@ -1354,6 +1354,8 @@ 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;
+    double deltaQ = length/static_cast<double>(numSteps-1);
     bool allDone = false;
     while (!allDone)
     {
@@ -1415,6 +1417,14 @@ namespace MDEvents
       // Is this box fully contained?
       if (verticesContained[i] >= maxVertices)
       {
+        coord_t boxCenter[nd];
+        box->getCenter(boxCenter);
+        // transform center of box to cylinder
+        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);
+        if (xchannel >= 0 || xchannel < numSteps ) signal_fit[xchannel] += box->getSignal();
         // Use the integrated sum of signal in the box
         signal += box->getSignal();
         errorSquared += box->getErrorSquared();
@@ -1437,7 +1447,7 @@ namespace MDEvents
         coord_t out[nd];
         radiusTransform.apply(boxCenter, out);
 
-        if (out[0] < diagonalSquared*0.72 + radius && std::fabs(out[1]) < diagonalSquared*0.72 + 0.5*length)
+        if (out[0] < std::sqrt(diagonalSquared*0.72 + radius*radius) && std::fabs(out[1]) < std::sqrt(diagonalSquared*0.72 + 0.25*length*length))
         {
           // If the center is closer than the size of the box, then it MIGHT be touching.
           // (We multiply by 0.72 (about sqrt(2)) to look for half the diagonal).