diff --git a/Framework/Algorithms/src/SofQW.cpp b/Framework/Algorithms/src/SofQW.cpp index 73d21353bd9a8a416ea21a86e15798e3ef403747..2d80eeacebdd0a9f0a088d878b023f745d445536 100644 --- a/Framework/Algorithms/src/SofQW.cpp +++ b/Framework/Algorithms/src/SofQW.cpp @@ -146,6 +146,7 @@ void SofQW::exec() { * @param[in] binParams The bin parameters from the user * @param[out] newAxis The 'vertical' axis defined by the given * parameters + * @param[out] ebinParams The 'horizontal' axis parameters (optional) * @return A pointer to the newly-created workspace */ API::MatrixWorkspace_sptr diff --git a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp index 8ce4215b4542c206b5d78f64e6ba5a7a0e339608..351de53c504c55966a19d5967a7e60a8c224782f 100644 --- a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp +++ b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp @@ -402,6 +402,7 @@ void SofQWNormalisedPolygon::initAngularCachesPSD( * @param[in] binParams The bin parameters from the user * @param[out] newAxis The 'vertical' axis defined by the given * parameters + * @param[out] ebinParams The 'horizontal' axis parameters (optional) * @return A pointer to the newly-created workspace */ RebinnedOutput_sptr SofQWNormalisedPolygon::setUpOutputWorkspace( diff --git a/Framework/DataObjects/src/FractionalRebinning.cpp b/Framework/DataObjects/src/FractionalRebinning.cpp index 0f45a52901e3d1b90421b4e9d572f0ff071110a0..f4ed9c474b619b26fdb9fb6aa541c8df631d3744 100644 --- a/Framework/DataObjects/src/FractionalRebinning.cpp +++ b/Framework/DataObjects/src/FractionalRebinning.cpp @@ -112,7 +112,7 @@ void calcRectangleIntersections( MatrixWorkspace_const_sptr outputWS, const std::vector<double> &yAxis, const Quadrilateral &inputQ, const size_t y_start, const size_t y_end, const size_t x_start, const size_t x_end, - std::vector<std::tuple<int, int, double>> &areaInfo) { + std::vector<std::tuple<size_t, size_t, double>> &areaInfo) { const auto &xAxis = outputWS->x(0); std::vector<double> width; for (size_t xi = x_start; xi < x_end; ++xi) { @@ -162,7 +162,7 @@ void calcTrapezoidYIntersections( MatrixWorkspace_const_sptr outputWS, const std::vector<double> &yAxis, const Quadrilateral &inputQ, const size_t y_start, const size_t y_end, const size_t x_start, const size_t x_end, - std::vector<std::tuple<int, int, double>> &areaInfo) { + std::vector<std::tuple<size_t, size_t, double>> &areaInfo) { // The algorithm proceeds as follows: // 1. Determine the left/right bin boundaries on the x- (horizontal)-grid. // 2. Loop along x, for each 1-output-bin wide strip construct a new input Q. @@ -203,7 +203,6 @@ void calcTrapezoidYIntersections( } // Step 1 - construct the left/right bin lims on the lines of the y-grid. - // First get equations for the top and bottom diagonal lines const double NaN = std::numeric_limits<double>::quiet_NaN(); const double DBL_EPS = std::numeric_limits<double>::epsilon(); std::vector<double> leftLim(nx * (ny + 1), NaN); @@ -214,10 +213,12 @@ void calcTrapezoidYIntersections( auto y1_it = yAxis.begin() + y_end + 1; const size_t ymax = (y_end == yAxis.size()) ? ny : (ny + 1); if ((mTop >= 0 && mBot >= 0) || (mTop < 0 && mBot < 0)) { - // Shape: /=/ or \=\ - top diag is left, bot is right or vice versa. + // Diagonals in same direction: For a given x-parallel line, + // Left limit given by one diagonal, right limit given by other double left_val, right_val; for (size_t yj = 0; yj < ymax; ++yj) { const size_t yjx = yj * nx; + // First, find the far left/right limits, given by the inputQ if (mTop >= 0) { left_val = (yAxis[yj + y_start] - cTop) / mTop; right_val = (yAxis[yj + y_start] - cBot) / mBot; @@ -245,6 +246,7 @@ void calcTrapezoidYIntersections( right_it = x1_it - 1; rightLim[nx - 1 + yjx] = lr_x; } + // Now populate the bin boundaries in between if (left_it < right_it && right_it != x1_it) { for (auto x_it = left_it; x_it != right_it; ++x_it) { leftLim[li + 1 + yjx] = *x_it; @@ -253,7 +255,7 @@ void calcTrapezoidYIntersections( } } } else { - // Shape: <=| or |=> - top diag is upper l/r, bot is lower l/r + // In this case, the diagonals are all on one side or the other. const size_t y2 = std::upper_bound(y0_it, y1_it, (mTop >= 0) ? ll_y : lr_y) - y0_it; const size_t y3 = @@ -302,7 +304,7 @@ void calcTrapezoidYIntersections( size_t vertBits = 0; size_t yj0, yj1; for (size_t xi = x_start; xi < x_end; ++xi) { - size_t xj = xi - x_start; + const size_t xj = xi - x_start; // Define new 1-bin wide input quadrilateral if (xi > x_start) { nll = nlr; @@ -450,7 +452,7 @@ void calcGeneralIntersections( MatrixWorkspace_const_sptr outputWS, const std::vector<double> &yAxis, const Quadrilateral &inputQ, const size_t qstart, const size_t qend, const size_t x_start, const size_t x_end, - std::vector<std::tuple<int, int, double>> &areaInfo) { + std::vector<std::tuple<size_t, size_t, double>> &areaInfo) { const auto &xAxis = outputWS->x(0); ConvexPolygon intersectOverlap; for (size_t yi = qstart; yi < qend; ++yi) { @@ -606,7 +608,7 @@ void rebinToFractionalOutput(const Quadrilateral &inputQ, // defined as rectangular. If the inputQ is is also rectangular or // trapezoidal, a simpler/faster way of calculating the intersection area // of all or some bins can be used. - std::vector<std::tuple<int, int, double>> areaInfo; + std::vector<std::tuple<size_t, size_t, double>> areaInfo; const double inputQArea = inputQ.area(); const QuadrilateralType inputQType = getQuadrilateralType(inputQ); if (inputQType == Rectangle) { @@ -621,8 +623,8 @@ void rebinToFractionalOutput(const Quadrilateral &inputQ, } for (auto ai : areaInfo) { - const int xi = std::get<0>(ai); - const int yi = std::get<1>(ai); + const size_t xi = std::get<0>(ai); + const size_t yi = std::get<1>(ai); const double weight = std::get<2>(ai) / inputQArea; PARALLEL_CRITICAL(overlap) { outputWS->mutableY(yi)[xi] += signal * weight; diff --git a/docs/source/algorithms/SofQWNormalisedPolygon-v1.rst b/docs/source/algorithms/SofQWNormalisedPolygon-v1.rst index c4bbb3ff6638f68e9e6a75f468570a572381b255..2ff556913a5b67d6010ca3f94d71516a4f9c437e 100644 --- a/docs/source/algorithms/SofQWNormalisedPolygon-v1.rst +++ b/docs/source/algorithms/SofQWNormalisedPolygon-v1.rst @@ -12,25 +12,49 @@ Description Converts a 2D workspace from :ref:`units <Unit Factory>` of spectrum number/**energy transfer** to the intensity as a function of **momentum transfer** :math:`Q` -and **energy transfer** :math:`\Delta E`. The rebinning is done as a -weighted sum of overlapping polygons with -fractional area tracking. The result is stored in a new workspace type: -**RebinnedOutput**. The new workspace presents the data as the -fractional counts divided by the fractional area. The biggest -consequence of this method is that in places where there are no counts -and no acceptance (no fractional areas), **nan**\ -s will result. - +and **energy transfer** :math:`\Delta E`. + +.. figure:: /images/RebinnedOutput.png + :align: center + +As shown in the figure, the input grid (pink-shaded parallelopiped, +aligned in scattering angle and energy transfer) is not parallel to the +output grid (square grid, aligned in :math:`Q` and energy). This means +that the output bins will only ever partially overlap the input data. To +account for this, the signal :math:`Y` and errors :math:`E` in the output +bin is calculated as the sum of all the input bins which overlap the +output bin, weighted by the fractional overlap area :math:`F_i`: + +.. math:: Y^{\mathrm{out}} = (\sum_i Y^{\mathrm{in}}_i F_i) / \sum_i F_i +.. math:: E^{\mathrm{out}} = \sqrt{\sum_i (E^{\mathrm{in}}_i F_i)^2} / \sum_i F_i + +.. warning:: Note that because the output bins contain fractions of multiple + input bins, the errors calculated for each output bins are no longer + independent, and so cannot be combined in quadrature. This means that + rebinning, summing, or integrating the output of this algorithm will + give *incorrect error values* because those Mantid algorithms use the + quadrature formular and assume independent bins. The *signal*, on the + other hand should still be correct on rebinning. Unary operations, such + as scaling the signal will not encounter this problem. + The algorithm operates in non-PSD mode by default. This means that all azimuthal angles and widths are forced to zero. PSD mode will determine the azimuthal angles and widths from the instrument geometry. This mode -is activated by placing the following named parameter in a Parameter +is activated by placing the following named parameter in the instrument definition file: *detector-neighbour-offset*. The integer value of this parameter should be the number of pixels that separates two pixels at the same -vertical position in adjacent tubes. - - -See :ref:`algm-SofQWCentre` for centre-point binning or :ref:`algm-SofQWPolygon` -for simpler and less precise but faster binning strategies. +vertical position in adjacent tubes. Note that in both non-PSD and PSD +modes, the scattering angle widths are determined from the detector +geometry and may vary from detector to detector as defined by the +instrument definition files. + +See :ref:`algm-SofQWCentre` for centre-point binning or :ref:`algm-SofQWPolygon` +for simpler and less precise but faster binning strategies. The speed-up +is from ignoring the azimuthal positions of the detectors (as for the non-PSD +mode in this algorithm) but in addition, :ref:`algm-SofQWPolygon` treats +all detectors as being the same, and characterised by a single width in +scattering angle. Thereafter, it weights the signal and error by the fractional +overlap, but does not then scale the weighted sum by :math:`\sum_i F_i`. Usage ----- diff --git a/docs/source/algorithms/SofQWPolygon-v1.rst b/docs/source/algorithms/SofQWPolygon-v1.rst index 17c4c3c0327dc1715ab315fec32006f97ec241bb..8028ef3cc859063671848def542efe96a2e6706d 100644 --- a/docs/source/algorithms/SofQWPolygon-v1.rst +++ b/docs/source/algorithms/SofQWPolygon-v1.rst @@ -14,9 +14,30 @@ of spectrum number/**energy transfer** to the intensity as a function of momentum transfer :math:`Q` and energy transfer :math:`\Delta E`. -The rebinning is done as a weighted sum of overlapping polygons. See -:ref:`algm-SofQWCentre` for centre-point binning or :ref:`algm-SofQWNormalisedPolygon` for -more complex and precise (but slower) binning strategy. +The rebinning is done as a weighted sum of overlapping polygons. +The polygon in :math:`Q-\Delta E` space is calculated from the +energy bin boundaries and the detector scattering angle :math:`2\theta`. +The detectors (pixels) are assumed to be uniform, and characterised +by a single angular width :math:`\Delta2\theta`. The signal and error +of the rebinned data (in :math:`Q-\Delta E` space) is then the +sum of the contributing pixels in each bin weighted by their fractional +overlap area: + +.. math:: Y^{\mathrm{out}} = (\sum_i Y^{\mathrm{in}}_i F_i) +.. math:: E^{\mathrm{out}} = \sqrt{\sum_i (E^{\mathrm{in}}_i F_i)^2} + +Unlike the more precise :ref:`algm-SofQWNormalisedPolygon` +algorithm, the final counts is not weighted by the sum of fractional +areas which means that the signal will be underestimated where the +bins are smaller, for example at the edges of detectors. +However, like the other algorithm, the output workspace has bins which +are no longer independent. Thus subsequent rebinning (or integration) +of the output of the algorithm will give incorrect error values. + +See :ref:`algm-SofQWCentre` for centre-point binning. +Alternatively, see :ref:`algm-SofQWNormalisedPolygon` for a +more complex and precise (but slower) binning strategy, where the actual +detector shape is calculated to obtain the input polygon. Usage ----- diff --git a/docs/source/images/RebinnedOutput.png b/docs/source/images/RebinnedOutput.png new file mode 100644 index 0000000000000000000000000000000000000000..39e1a7c688435941fc70623fb6faf712097f5c18 Binary files /dev/null and b/docs/source/images/RebinnedOutput.png differ diff --git a/docs/source/release/v3.12.0/spectroscopy.rst b/docs/source/release/v3.12.0/spectroscopy.rst index 7839e668f7836233ec949792c25f4b5cdbbb69b9..c327a21767bbd88294b14be30836de85e2202f16 100644 --- a/docs/source/release/v3.12.0/spectroscopy.rst +++ b/docs/source/release/v3.12.0/spectroscopy.rst @@ -9,6 +9,9 @@ Spectroscopy Changes putting new features at the top of the section, followed by improvements, followed by bug fixes. +- The algorithms :ref:`algm-SofQWCentre`, :ref:`algm-SofQWPolygon` and :ref:`algm-SofQWNormalisedPolygon`, which rebin an inelastic workspace (has a `DeltaE` axis) from spectrum numbers (angle) to `MomentumTransfer` may now rebin the energy (`DeltaE`) axis as well as the :math:`|Q|` (`MomentumTransfer`) axes. +- :ref:`algm-SofQWNormalisedPolygon` now has uses a faster method for calculating the polygon intersections. + Direct Geometry ---------------