MDNormalization.cpp

      double fmom = (kfmax - kfmin) / (kEnd - kStart);
      double fh = (hEnd - hStart) / (kEnd - kStart);
      double fl = (lEnd - lStart) / (kEnd - kStart);
      for (size_t i = 0; i < kNBins; i++) {
        double ki = m_kX[i];
        if (((kStart - ki) * (kEnd - ki) < 0)) {
          // if ki is between kStart and kEnd, then hi and li will be between
          // hStart, hEnd and lStart, lEnd and momi will be between kfmin and
          // kfmax
          double hi = fh * (ki - kStart) + hStart;
          double li = fl * (ki - kStart) + lStart;
          if ((hi >= m_hX[0]) && (hi <= m_hX[hNBins-1]) && (li >= m_lX[0]) &&
              (li <= m_lX[lNBins-1])) {
            double momi = fmom * (ki - kStart) + kfmin;
            intersections.push_back({{hi, ki, li, momi}});
          }
        }
      }
    }

    // calculate intersections with planes perpendicular to l
    if (fabs(lStart - lEnd) > eps) {
      double fmom = (kfmax - kfmin) / (lEnd - lStart);
      double fh = (hEnd - hStart) / (lEnd - lStart);
      double fk = (kEnd - kStart) / (lEnd - lStart);

      for (size_t i = 0; i < lNBins; i++) {
        double li = m_lX[i];
        if (((lStart - li) * (lEnd - li) < 0)) {
          double hi = fh * (li - lStart) + hStart;
          double ki = fk * (li - lStart) + kStart;
          if ((hi >= m_hX[0]) && (hi <= m_hX[hNBins-1]) && (ki >= m_kX[0]) &&
               (ki <= m_kX[kNBins-1])) {
            double momi = fmom * (li - lStart) + kfmin;
            intersections.push_back({{hi, ki, li, momi}});
          }
        }
      }
    }
    // intersections with dE
    if (!m_dEIntegrated) {
      for (size_t i = 0; i < eNBins; i++) {
        double kfi = m_eX[i];
        if ((kfi - kfmin) * (kfi - kfmax) <= 0) {
          double h = qin.X() - qout.X() * kfi;
          double k = qin.Y() - qout.Y() * kfi;
          double l = qin.Z() - qout.Z() * kfi;
          if ((h >= m_hX[0]) && (h <= m_hX[hNBins-1]) && (k >= m_kX[0]) && (k <= m_kX[kNBins-1]) &&
              (l >= m_lX[0]) && (l <= m_lX[lNBins-1])) {
            intersections.push_back({{h, k, l, kfi}});
          }
        }
      }
    }

    // endpoints
    if ((hStart >= m_hX[0]) && (hStart <= m_hX[hNBins-1]) && (kStart >= m_kX[0]) &&
        (kStart <= m_kX[kNBins-1]) && (lStart >= m_lX[0]) && (lStart <= m_lX[lNBins-1])) {
      intersections.push_back({{hStart, kStart, lStart, kfmin}});
    }
    if ((hEnd >= m_hX[0]) && (hEnd <= m_hX[hNBins-1]) && (kEnd >= m_kX[0]) &&
        (kEnd <= m_kX[kNBins-1]) && (lEnd >= m_lX[0]) && (lEnd <= m_lX[lNBins-1])) {
      intersections.push_back({{hEnd, kEnd, lEnd, kfmax}});
    }

    // sort intersections by final momentum
    std::stable_sort(intersections.begin(), intersections.end(), compareMomentum);
}

void MDNormalization::calcIntegralsForIntersections(const std::vector<double> &xValues, const API::MatrixWorkspace &integrFlux, size_t sp, std::vector<double> &yValues)
{
    assert(xValues.size() == yValues.size());

    // the x-data from the workspace
    const auto &xData = integrFlux.x(sp);
    const double xStart = xData.front();
    const double xEnd = xData.back();

    // the values in integrFlux are expected to be integrals of a non-negative
    // function
    // ie they must make a non-decreasing function
    const auto &yData = integrFlux.y(sp);
    size_t spSize = yData.size();

    const double yMin = 0.0;
    const double yMax = yData.back();

    size_t nData = xValues.size();
    // all integrals below xStart must be 0
    if (xValues[nData - 1] < xStart) {
      std::fill(yValues.begin(), yValues.end(), yMin);
      return;
    }

    // all integrals above xEnd must be equal tp yMax
    if (xValues[0] > xEnd) {
      std::fill(yValues.begin(), yValues.end(), yMax);
      return;
    }

    size_t i = 0;
    // integrals below xStart must be 0
    while (i < nData - 1 && xValues[i] < xStart) {
      yValues[i] = yMin;
      i++;
    }
    size_t j = 0;
    for (; i < nData; i++) {
      // integrals above xEnd must be equal tp yMax
      if (j >= spSize - 1) {
        yValues[i] = yMax;
      } else {
        double xi = xValues[i];
        while (j < spSize - 1 && xi > xData[j])
          j++;
        // if x falls onto an interpolation point return the corresponding y
        if (xi == xData[j]) {
          yValues[i] = yData[j];
        } else if (j == spSize - 1) {
          // if we get above xEnd it's yMax
          yValues[i] = yMax;
        } else if (j > 0) {
          // interpolate between the consecutive points
          double x0 = xData[j - 1];
          double x1 = xData[j];
          double y0 = yData[j - 1];
          double y1 = yData[j];
          yValues[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
        } else // j == 0
        {
          yValues[i] = yMin;
        }
      }
    }
}

} // namespace MDAlgorithms
} // namespace Mantid