MDHistoWorkspace.cpp

    signal_t a = m_signals[i];
    signal_t da2 = m_errorsSquared[i];
    if (a <= 0) {
      m_signals[i] = filler;
      m_errorsSquared[i] = 0;
    } else {
      m_signals[i] = std::log10(a);
      m_errorsSquared[i] = 0.1886117 * da2 / (a * a); // 0.1886117  = ln(10)^-2
    }
  }
}

//----------------------------------------------------------------------------------------------
/** Perform the exp() function on each signal in the workspace.
 *
 * Error propagation of \f$ f = exp(a) \f$  is given by:
 * \f$ df^2 = f^2 * da^2 \f$
 */
void MDHistoWorkspace::exp() {
  for (size_t i = 0; i < m_length; ++i) {
    signal_t f = std::exp(m_signals[i]);
    signal_t da2 = m_errorsSquared[i];
    m_signals[i] = f;
    m_errorsSquared[i] = f * f * da2;
  }
}

//----------------------------------------------------------------------------------------------
/** Perform the power function (signal^exponent) on each signal S in the
 *workspace.
 *
 * Error propagation of \f$ f = a^b \f$  is given by:
 * \f$ df^2 = f^2 * b^2 * (da^2 / a^2) \f$
 */
void MDHistoWorkspace::power(double exponent) {
  double exponent_squared = exponent * exponent;
  for (size_t i = 0; i < m_length; ++i) {
    signal_t a = m_signals[i];
    signal_t f = std::pow(a, exponent);
    signal_t da2 = m_errorsSquared[i];
    m_signals[i] = f;
    m_errorsSquared[i] = f * f * exponent_squared * da2 / (a * a);
  }
}

//==============================================================================================
//============================== BOOLEAN OPERATIONS
//============================================
//==============================================================================================

//----------------------------------------------------------------------------------------------
/** A boolean &= (and) operation, element-by-element, for two
 *MDHistoWorkspace's.
 *
 * 0.0 is "false", all other values are "true". All errors are set to 0.
 *
 * @param b :: workspace on the RHS of the operation
 * @return *this after operation */
MDHistoWorkspace &MDHistoWorkspace::operator&=(const MDHistoWorkspace &b) {
  checkWorkspaceSize(b, "&= (and)");
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = ((m_signals[i] != 0 && !m_masks[i]) &&
                    (b.m_signals[i] != 0 && !b.m_masks[i]))
                       ? 1.0
                       : 0.0;
    m_errorsSquared[i] = 0;
  }
  return *this;
}

//----------------------------------------------------------------------------------------------
/** A boolean |= (or) operation, element-by-element, for two MDHistoWorkspace's.
 *
 * 0.0 is "false", all other values are "true". All errors are set to 0.
 *
 * @param b :: workspace on the RHS of the operation
 * @return *this after operation */
MDHistoWorkspace &MDHistoWorkspace::operator|=(const MDHistoWorkspace &b) {
  checkWorkspaceSize(b, "|= (or)");
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = ((m_signals[i] != 0 && !m_masks[i]) ||
                    (b.m_signals[i] != 0 && !b.m_masks[i]))
                       ? 1.0
                       : 0.0;
    m_errorsSquared[i] = 0;
  }
  return *this;
}

//----------------------------------------------------------------------------------------------
/** A boolean ^= (xor) operation, element-by-element, for two
 *MDHistoWorkspace's.
 *
 * 0.0 is "false", all other values are "true". All errors are set to 0.
 *
 * @param b :: workspace on the RHS of the operation
 * @return *this after operation */
MDHistoWorkspace &MDHistoWorkspace::operator^=(const MDHistoWorkspace &b) {
  checkWorkspaceSize(b, "^= (xor)");
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = ((m_signals[i] != 0 && !m_masks[i]) ^
                    (b.m_signals[i] != 0 && !b.m_masks[i]))
                       ? 1.0
                       : 0.0;
    m_errorsSquared[i] = 0;
  }
  return *this;
}

//----------------------------------------------------------------------------------------------
/** A boolean not operation, performed in-place.
 * All errors are set to 0.
 *
 * 0.0 is "false", all other values are "true". All errors are set to 0.
 */
void MDHistoWorkspace::operatorNot() {
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = (m_signals[i] == 0.0 || m_masks[i]);
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is < b[i].
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param b :: workspace on the RHS of the comparison.
 */
void MDHistoWorkspace::lessThan(const MDHistoWorkspace &b) {
  checkWorkspaceSize(b, "lessThan");
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = (m_signals[i] < b.m_signals[i]) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is < signal.
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param signal :: signal value on the RHS of the comparison.
 */
void MDHistoWorkspace::lessThan(const signal_t signal) {
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = (m_signals[i] < signal) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is > b[i].
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param b :: workspace on the RHS of the comparison.
 */
void MDHistoWorkspace::greaterThan(const MDHistoWorkspace &b) {
  checkWorkspaceSize(b, "greaterThan");
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = (m_signals[i] > b.m_signals[i]) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is > signal.
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param signal :: signal value on the RHS of the comparison.
 */
void MDHistoWorkspace::greaterThan(const signal_t signal) {
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = (m_signals[i] > signal) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is == b[i].
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param b :: workspace on the RHS of the comparison.
 * @param tolerance :: accept this deviation from a perfect equality
 */
void MDHistoWorkspace::equalTo(const MDHistoWorkspace &b,
                               const signal_t tolerance) {
  checkWorkspaceSize(b, "equalTo");
  for (size_t i = 0; i < m_length; ++i) {
    signal_t diff = fabs(m_signals[i] - b.m_signals[i]);
    m_signals[i] = (diff < tolerance) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Turn this workspace into a boolean workspace, where
 * signal[i] -> becomes true (1.0) if it is == signal.
 * signal[i] -> becomes false (0.0) otherwise
 * Errors are set to 0.
 *
 * @param signal :: signal value on the RHS of the comparison.
 * @param tolerance :: accept this deviation from a perfect equality
 */
void MDHistoWorkspace::equalTo(const signal_t signal,
                               const signal_t tolerance) {
  for (size_t i = 0; i < m_length; ++i) {
    signal_t diff = fabs(m_signals[i] - signal);
    m_signals[i] = (diff < tolerance) ? 1.0 : 0.0;
    m_errorsSquared[i] = 0;
  }
}

//----------------------------------------------------------------------------------------------
/** Copy the values from another workspace onto this workspace, but only
 * where a mask is true (non-zero)
 *
 * For example, in matlab or numpy python, you might write something like:
 *  "mask = (array < 5.0); array[mask] = other[mask];"
 *
 * The equivalent here is:
 *  mask = array;
 *  mask.lessThan(5.0);
 *  array.setUsingMask(mask, other);
 *
 * @param mask :: MDHistoWorkspace where (signal == 0.0) means false, and
 *(signal != 0.0) means true.
 * @param values :: MDHistoWorkspace of values to copy.
 */
void MDHistoWorkspace::setUsingMask(const MDHistoWorkspace &mask,
                                    const MDHistoWorkspace &values) {
  checkWorkspaceSize(mask, "setUsingMask");
  checkWorkspaceSize(values, "setUsingMask");
  for (size_t i = 0; i < m_length; ++i) {
    if (mask.m_signals[i] != 0.0) {
      m_signals[i] = values.m_signals[i];
      m_errorsSquared[i] = values.m_errorsSquared[i];
    }
  }
}

//----------------------------------------------------------------------------------------------
/** Copy the values from another workspace onto this workspace, but only
 * where a mask is true (non-zero)
 *
 * For example, in matlab or numpy python, you might write something like:
 *  "mask = (array < 5.0); array[mask] = other[mask];"
 *
 * The equivalent here is:
 *  mask = array;
 *  mask.lessThan(5.0);
 *  array.setUsingMask(mask, other);
 *
 * @param mask :: MDHistoWorkspace where (signal == 0.0) means false, and
 *(signal != 0.0) means true.
 * @param signal :: signal to set everywhere mask is true
 * @param error :: error (not squared) to set everywhere mask is true
 */
void MDHistoWorkspace::setUsingMask(const MDHistoWorkspace &mask,
                                    const signal_t signal,
                                    const signal_t error) {
  signal_t errorSquared = error * error;
  checkWorkspaceSize(mask, "setUsingMask");
  for (size_t i = 0; i < m_length; ++i) {
    if (mask.m_signals[i] != 0.0) {
      m_signals[i] = signal;
      m_errorsSquared[i] = errorSquared;
    }
  }
}

/**
Setter for the masking region.
Does not perform any clearing. Multiple calls are compounded.
@param maskingRegion : Implicit function defining mask region.
*/
void MDHistoWorkspace::setMDMasking(
    Mantid::Geometry::MDImplicitFunction *maskingRegion) {
  if (maskingRegion != nullptr) {
    for (size_t i = 0; i < this->getNPoints(); ++i) {
      // If the function masks the point, then mask it, otherwise leave it as it
      // is.
      if (maskingRegion->isPointContained(this->getCenter(i))) {
        this->setMDMaskAt(i, true);
      }
    }
    delete maskingRegion;
  }
}

/**
 * Set the masking
 * @param index : linear index to mask
 * @param mask : True to mask. False to clear.
 */
void MDHistoWorkspace::setMDMaskAt(const size_t &index, bool mask) {
  m_masks[index] = mask;
  if (mask) {
    // Set signal and error of masked points to the value of MDMaskValue
    this->setSignalAt(index, MDMaskValue);
    this->setErrorSquaredAt(index, MDMaskValue);
  }
}

/**
 * Clear any existing masking.
 * Note that this clears the mask flag but does not restore the data
 * which was set to NaN when it was masked.
 */
void MDHistoWorkspace::clearMDMasking() {
  for (size_t i = 0; i < this->getNPoints(); ++i) {
    m_masks[i] = false;
  }
}

uint64_t MDHistoWorkspace::getNEvents() const {
  volatile uint64_t cach = this->m_nEventsContributed;
  if (cach != this->m_nEventsContributed) {
    if (!m_numEvents)
      m_nEventsContributed = std::numeric_limits<uint64_t>::quiet_NaN();
    else
      m_nEventsContributed = sumNContribEvents();
  }
  return m_nEventsContributed;
}

uint64_t MDHistoWorkspace::sumNContribEvents() const {
  uint64_t sum(0);
  for (size_t i = 0; i < m_length; ++i)
    sum += uint64_t(m_numEvents[i]);

  return sum;
}

/**
 * Get the Q frame system (if any) to use.
*/
// clang-format off
GCC_DIAG_OFF(strict-aliasing)
// clang-format on
Kernel::SpecialCoordinateSystem
MDHistoWorkspace::getSpecialCoordinateSystem() const {
  MDFramesToSpecialCoordinateSystem converter;
  auto coordinatesFromMDFrames = converter(this);
  auto coordinates = m_coordSystem;
  if (coordinatesFromMDFrames) {
    coordinates = coordinatesFromMDFrames.get();
  }
  return coordinates;
}

/**
Set the special coordinate system (if any) to use.
@param coordinateSystem : Special coordinate system to use.
*/
void MDHistoWorkspace::setCoordinateSystem(
    const Kernel::SpecialCoordinateSystem coordinateSystem) {
  m_coordSystem = coordinateSystem;
}

/**
 * Static helper method.
 * @return The size of an element in the MDEventWorkspace.
 */
size_t MDHistoWorkspace::sizeOfElement() {
  return (3 * sizeof(signal_t)) + sizeof(bool);
}

/**
Preferred normalization to use for visual purposes.
*/
MDNormalization MDHistoWorkspace::displayNormalization() const {
  return m_displayNormalization; // Normalize by the number of events.
}

/**
Preferred normalization to use for visual purposes.
*/
MDNormalization MDHistoWorkspace::displayNormalizationHisto() const {
  return displayNormalization(); // Normalize by the number of events.
}

void MDHistoWorkspace::setDisplayNormalization(
    const Mantid::API::MDNormalization &preferredNormalization) {
  m_displayNormalization = preferredNormalization;
}

} // namespace Mantid
} // namespace DataObjects