Refs #4202. Fix a factor of 2 error in bin size calc.

mantidproject · Jan 26, 2012 · 93b74bf · 93b74bf
1 parent a7c7232
commit 93b74bf
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 17 deletions.
diff --git a/Code/Mantid/Framework/Algorithms/src/SofQW2.cpp b/Code/Mantid/Framework/Algorithms/src/SofQW2.cpp
@@ -87,13 +87,15 @@ namespace Mantid
         qCalculator = &SofQW2::calculateIndirectQ;
       }
 
-
+      PARALLEL_FOR2(inputWS, outputWS)
       for(int64_t i = 0; i < static_cast<int64_t>(nTheta); ++i) // signed for openmp
       {
+        PARALLEL_START_INTERUPT_REGION
+
         const double theta = m_thetaPts[i];
         const double efixed = getEFixed(inputWS->getDetector(i));
-        const double thetaLower = theta - m_thetaWidth;
-        const double thetaUpper = theta + m_thetaWidth;
+        const double thetaLower = theta - 0.5*m_thetaWidth;
+        const double thetaUpper = theta + 0.5*m_thetaWidth;
 
         for(size_t j = 0; j < nenergyBins; ++j)
         {
@@ -111,7 +113,24 @@ namespace Mantid
 
           rebinToOutput(inputQ, inputWS, i, j, outputWS);
         }
+
+        PARALLEL_END_INTERUPT_REGION
+      }
+      PARALLEL_CHECK_INTERUPT_REGION
+
+      const size_t nOutputHist(outputWS->getNumberHistograms());
+      // The errors need square-rooting
+      PARALLEL_FOR1(outputWS)
+      for(int64_t i = 0; i < static_cast<int64_t>(nOutputHist); ++i)
+      {
+        PARALLEL_START_INTERUPT_REGION
+
+        MantidVec& errors = outputWS->dataE(i);
+        std::transform(errors.begin(), errors.end(), errors.begin(), (double (*)(double))std::sqrt);
+
+        PARALLEL_END_INTERUPT_REGION
       }
+      PARALLEL_CHECK_INTERUPT_REGION
 
       // Divide by the bin width if the input is a distribution
       if( inputWS->isDistribution() )
@@ -146,7 +165,9 @@ namespace Mantid
           try
           {
             ConvexPolygon overlap = intersectionByLaszlo(outputQ, inputQ);
-            outputWS->dataY(qi)[ei] += inputWS->readY(i)[j] * (overlap.area()/inputQ.area());
+            const double weight = overlap.area()/inputQ.area();
+            outputWS->dataY(qi)[ei] += inputWS->readY(i)[j] * weight;
+            outputWS->dataE(qi)[ei] += std::pow(inputWS->readE(i)[j] * weight, 2);
           }
           catch(Geometry::NoIntersectionException &)
           {}

diff --git a/Code/Mantid/Framework/Algorithms/test/SofQW2Test.h b/Code/Mantid/Framework/Algorithms/test/SofQW2Test.h
@@ -69,28 +69,28 @@ class SofQW2Test : public CxxTest::TestSuite
 //    std::cerr << "\n";
 //    for(int i = 0; i < 6; ++i)
 //    {
-//      std::cerr << std::setprecision(12) << result->readY(i)[bins[i]] << "  " << result->readE(i)[bins[i]];
+//      std::cerr << std::setprecision(12) << result->readY(i)[bins[i]] << "  " << result->readE(i)[bins[i]] << "\n";
 //    }
 //    std::cerr << "\n";
 
     const double delta(1e-08);
-    TS_ASSERT_DELTA( result->readY(0)[1160], 5.41214271, delta);
-    //TS_ASSERT_DELTA( result->readE(0)[1160], 0.062926550, delta);
+    TS_ASSERT_DELTA( result->readY(0)[1160], 0.0, delta);
+    TS_ASSERT_DELTA( result->readE(0)[1160], 0.0, delta);
 
-    TS_ASSERT_DELTA( result->readY(1)[1145], 13.12828673, delta);
-    //TS_ASSERT_DELTA( result->readE(1)[1145], 0.040755760, delta);
+    TS_ASSERT_DELTA( result->readY(1)[1145], 14.21487562, delta);
+    TS_ASSERT_DELTA( result->readE(1)[1145], 0.27891342, delta);
 
-    TS_ASSERT_DELTA( result->readY(2)[1200], 6.63459775, delta);
-    //TS_ASSERT_DELTA( result->readE(2)[1200], 0.021459611, delta);
+    TS_ASSERT_DELTA( result->readY(2)[1200], 5.33576534, delta);
+    TS_ASSERT_DELTA( result->readE(2)[1200], 0.16793570, delta);
 
-    TS_ASSERT_DELTA( result->readY(3)[99], 0.17718250  , delta);
-    //TS_ASSERT_DELTA( result->readE(3)[99],   0.007609629, delta);
+    TS_ASSERT_DELTA( result->readY(3)[99],  0.17843416 , delta);
+    TS_ASSERT_DELTA( result->readE(3)[99],  0.07401977, delta);
 
-    TS_ASSERT_DELTA( result->readY(4)[1654],0.07161639 , delta);
-    //TS_ASSERT_DELTA( result->readE(4)[1654], 0.00182863, delta);
+    TS_ASSERT_DELTA( result->readY(4)[1654], 0.06845460, delta);
+    TS_ASSERT_DELTA( result->readE(4)[1654], 0.02002920, delta);
 
-    TS_ASSERT_DELTA( result->readY(5)[1025], 0.14990359, delta);
-    //TS_ASSERT_DELTA( result->readE(5)[1025], 0.00410309, delta);
+    TS_ASSERT_DELTA( result->readY(5)[1025], 0.20644528, delta);
+    TS_ASSERT_DELTA( result->readE(5)[1025], 0.04115725, delta);
 
     AnalysisDataService::Instance().remove(inputWS);
     AnalysisDataService::Instance().remove(outputWS);