diff --git a/Framework/MDAlgorithms/src/Integrate3DEvents.cpp b/Framework/MDAlgorithms/src/Integrate3DEvents.cpp index f8c9dcd723f92093bed26c896d08ec78062c738d..caa82bda406d39eecd555c2ff8ed23e9af8c6aa7 100644 --- a/Framework/MDAlgorithms/src/Integrate3DEvents.cpp +++ b/Framework/MDAlgorithms/src/Integrate3DEvents.cpp @@ -181,9 +181,9 @@ Integrate3DEvents::integrateWeakPeak( const auto &events = *result; const auto &directions = shape->directions(); - const auto &abcBackgroundInnerRadii = shape->abcRadiiBackgroundInner(); - const auto &abcBackgroundOuterRadii = shape->abcRadiiBackgroundOuter(); - const auto &abcRadii = shape->abcRadii(); + auto abcBackgroundInnerRadii = shape->abcRadiiBackgroundInner(); + auto abcBackgroundOuterRadii = shape->abcRadiiBackgroundOuter(); + auto abcRadii = shape->abcRadii(); const auto max_sigma = std::get<2>(libPeak); auto rValues = calculateRadiusFactors(params, max_sigma); @@ -192,6 +192,8 @@ Integrate3DEvents::integrateWeakPeak( correctForDetectorEdges(rValues, params.E1Vectors, center, abcRadii, abcBackgroundInnerRadii, abcBackgroundOuterRadii); + + if (!isPeakOnDetector) return shape; @@ -209,7 +211,6 @@ Integrate3DEvents::integrateWeakPeak( const auto fracError = std::get<1>(libPeak); inti = peak_w_back - ratio * backgrd; - sigi = inti + ratio * ratio * backgrd; // correct for fractional intensity sigi = sigi / pow(inti, 2); @@ -218,7 +219,17 @@ Integrate3DEvents::integrateWeakPeak( inti = inti * frac; sigi = sqrt(sigi) * inti; - return shape; + // scale integration shape by fractional amount + for (size_t i = 0; i < abcRadii.size(); ++i) { + abcRadii[i] *= frac; + abcBackgroundInnerRadii[i] *= frac; + abcBackgroundOuterRadii[i] *= frac; + } + + return boost::make_shared<const PeakShapeEllipsoid>( + shape->directions(), abcRadii, abcBackgroundInnerRadii, + abcBackgroundOuterRadii, Mantid::Kernel::QLab, + "IntegrateEllipsoidsTwoStep"); } double Integrate3DEvents::estimateSignalToNoiseRatio( @@ -267,7 +278,7 @@ double Integrate3DEvents::estimateSignalToNoiseRatio( double peak_w_back = numInEllipsoid(events, eigen_vectors, peakRadii); double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3)); - double inti = peak_w_back - ratio * backgrd; + auto inti = peak_w_back - ratio * backgrd; auto sigi = sqrt(peak_w_back + ratio * ratio * backgrd); return inti / sigi; diff --git a/Framework/MDAlgorithms/test/Integrate3DEventsTest.h b/Framework/MDAlgorithms/test/Integrate3DEventsTest.h index 2d062d3ec01f0daa63df85854bd6bf3a1643a320..94132c67a8c9a9c17aa09c8af199c51f69b447c8 100644 --- a/Framework/MDAlgorithms/test/Integrate3DEventsTest.h +++ b/Framework/MDAlgorithms/test/Integrate3DEventsTest.h @@ -185,8 +185,8 @@ public: // to be weighted by the fraction of strong peak contained in a standard // core. This is not exactly the same because of the weighting from the // strong peak - TS_ASSERT_DELTA(weak_inti, 83.6960, 0.5); - TS_ASSERT_DELTA(weak_sigi, 8.37, 0.1); + TS_ASSERT_DELTA(weak_inti, numWeakEvents, 0.5); + TS_ASSERT_DELTA(weak_sigi, 10, 0.1); weak_inti = 0; weak_sigi = 0; @@ -260,7 +260,7 @@ public: // core. This is not exactly the same because of the weighting from the // strong peak TS_ASSERT_DELTA(weak_inti, numWeakEvents, 35); - TS_ASSERT_DELTA(weak_sigi, 8.62, 0.2); + TS_ASSERT_DELTA(weak_sigi, 10, 0.2); } void test_estimateSignalToNoiseRatioInPerfectCase() { @@ -298,17 +298,17 @@ public: const auto ratio2 = integrator.estimateSignalToNoiseRatio(params, peak_2); const auto ratio3 = integrator.estimateSignalToNoiseRatio(params, peak_3); - TS_ASSERT_DELTA(ratio1, numStrongEvents, 0.0001); - TS_ASSERT_DELTA(ratio2, numWeakEvents, 0.0001); - TS_ASSERT_DELTA(ratio3, numWeakEvents / 2, 0.0001); + TS_ASSERT_DELTA(ratio1, numStrongEvents / 100, 1e-4); + TS_ASSERT_DELTA(ratio2, numWeakEvents / 10, 1e-4); + TS_ASSERT_DELTA(ratio3, 7.071, 1e-4); } void test_estimateSignalToNoiseRatioWithBackgroundAndOnePercentCulling() { - doTestSignalToNoiseRatio(true, 171.90, 1.2632, 0.1824); + doTestSignalToNoiseRatio(true, 99.3898, 5.4788, 1.0597); } void test_estimateSignalToNoiseRatioWithBackgroundAndNoOnePercentCulling() { - doTestSignalToNoiseRatio(false, 160.33, 1.088, 0.094); + doTestSignalToNoiseRatio(false, 99.3417, 5.0972, 0.5821); } private: