C++ MatrixWorkspace_const_sptr::getNumberHistograms方法代码示例

    /**
     * Return true if the two workspaces are compatible for this operation
     * Virtual: will be overridden as needed.
     * @param lhs :: left-hand workspace to check
     * @param rhs :: right-hand workspace to check
     * @return flag for the compatibility to the two workspaces
     */
    bool BinaryOperation::checkCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const
    {
      Unit_const_sptr lhs_unit;
      Unit_const_sptr rhs_unit;
      if ( lhs->axes() && rhs->axes() ) // If one of these is a WorkspaceSingleValue then we don't want to check units match
      {
        lhs_unit = lhs->getAxis(0)->unit();
        rhs_unit = rhs->getAxis(0)->unit();
      }

      const std::string lhs_unitID = ( lhs_unit ? lhs_unit->unitID() : "" );
      const std::string rhs_unitID = ( rhs_unit ? rhs_unit->unitID() : "" );

      // Check the workspaces have the same units and distribution flag
      if ( lhs_unitID != rhs_unitID && lhs->blocksize() > 1 && rhs->blocksize() > 1 )
      {
        g_log.error("The two workspace are not compatible because they have different units on the X axis.");
        return false;
      }

      // Check the size compatibility
      if (!checkSizeCompatibility(lhs,rhs))
      {
        std::ostringstream ostr;
        ostr<<"The sizes of the two workspaces " <<
            "(" << lhs->getName() << ": " << lhs->getNumberHistograms() << " spectra, blocksize " << lhs->blocksize() << ")"
            << " and " <<
            "(" << rhs->getName() << ": " << rhs->getNumberHistograms() << " spectra, blocksize " << rhs->blocksize() << ")"
            << " are not compatible for algorithm "<<this->name();
        g_log.error() << ostr.str() << std::endl;
        throw std::invalid_argument( ostr.str() );
      }

      return true;
    }

/** Checks if workspaces input to Q1D or Qxy are reasonable
  @param dataWS data workspace
  @param binAdj (WavelengthAdj) workpace that will be checked to see if it has
  one spectrum and the same number of bins as dataWS
  @param detectAdj (PixelAdj) passing NULL for this wont raise an error, if set
  it will be checked this workspace has as many histograms as dataWS each with
  one bin
  @param qResolution: the QResolution workspace
  @throw invalid_argument if the workspaces are not mututially compatible
*/
void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
                           API::MatrixWorkspace_const_sptr binAdj,
                           API::MatrixWorkspace_const_sptr detectAdj,
                           API::MatrixWorkspace_const_sptr qResolution) {

  // Check the compatibility of dataWS, binAdj and detectAdj
  examineInput(dataWS, binAdj, detectAdj);

  // Check the compatibility of the QResolution workspace
  if (qResolution) {
    // We require the same number of histograms
    if (qResolution->getNumberHistograms() != dataWS->getNumberHistograms()) {
      throw std::invalid_argument("The QResolution should have one spectrum"
                                  "per spectrum of the input workspace");
    }

    // We require the same binning for the input workspace and the q resolution
    // workspace
    MantidVec::const_iterator reqX = dataWS->readX(0).begin();
    MantidVec::const_iterator qResX = qResolution->readX(0).begin();
    for (; reqX != dataWS->readX(0).end(); ++reqX, ++qResX) {
      if (*reqX != *qResX) {
        throw std::invalid_argument(
            "The QResolution needs to have the same binning as"
            "as the input workspace.");
      }
    }
  }
}

/** Performs a simple check to see if the sizes of two workspaces are compatible
 *for a binary operation
 *  In order to be size compatible then the larger workspace
 *  must divide be the size of the smaller workspace leaving no remainder
 *
 *  @param lhs :: the first workspace to compare
 *  @param rhs :: the second workspace to compare
 *  @retval true The two workspaces are size compatible
 *  @retval false The two workspaces are NOT size compatible
 */
std::string Multiply::checkSizeCompatibility(
    const API::MatrixWorkspace_const_sptr lhs,
    const API::MatrixWorkspace_const_sptr rhs) const {
  if (!m_keepEventWorkspace && !m_AllowDifferentNumberSpectra) {
    // Fallback on the default checks
    return CommutativeBinaryOperation::checkSizeCompatibility(lhs, rhs);
  } else {

    // A SingleValueWorkspace on the right, or matches anything
    if (rhs->size() == 1)
      return "";

    // A SingleValueWorkspace on the left only matches if rhs was single value
    // too. Why are you using mantid to do simple math?!?
    if (lhs->size() == 1)
      return "The left side cannot contain a single value if the right side "
             "isn't also a single value.";

    // RHS only has one value (1D vertical), so the number of histograms needs
    // to match.
    // Each lhs spectrum will be divided by that scalar
    if (rhs->blocksize() == 1 &&
        lhs->getNumberHistograms() == rhs->getNumberHistograms())
      return "";

    if (m_matchXSize) {
      // Past this point, for a 2D WS operation, we require the X arrays to
      // match. Note this only checks the first spectrum
      if (!WorkspaceHelpers::matchingBins(*lhs, *rhs, true)) {
        return "X arrays must match when multiplying 2D workspaces.";
      }
    }

    // We don't need to check for matching bins for events. Yay events!
    const size_t rhsSpec = rhs->getNumberHistograms();

    // If the rhs has a single spectrum, then we can divide. The block size does
    // NOT need to match,
    if (rhsSpec == 1)
      return "";

    // Are we allowing the division by different # of spectra, using detector
    // IDs to match up?
    if (m_AllowDifferentNumberSpectra) {
      return "";
    }

    // Otherwise, the number of histograms needs to match, but the block size of
    // each does NOT need to match.

    if (lhs->getNumberHistograms() == rhs->getNumberHistograms()) {
      return "";
    } else {
      return "Number of histograms not identical.";
    }
  }
}

/** Initialization method:
@param bkgWS    -- shared pointer to the workspace which contains background
@param sourceWS -- shared pointer to the workspace to remove background from
@param emode    -- energy conversion mode used during internal units conversion
(0 -- elastic, 1-direct, 2 indirect, as defined in Units conversion
@param pLog     -- pointer to the logger class which would report errors
@param nThreads -- number of threads to be used for background removal
@param inPlace  -- if the background removal occurs from the existing workspace
or target workspace has to be cloned.
*/
void BackgroundHelper::initialize(const API::MatrixWorkspace_const_sptr &bkgWS,
                                  const API::MatrixWorkspace_sptr &sourceWS,
                                  int emode, Kernel::Logger *pLog, int nThreads,
                                  bool inPlace) {
  m_bgWs = bkgWS;
  m_wkWS = sourceWS;
  m_Emode = emode;
  m_pgLog = pLog;
  m_inPlace = inPlace;

  std::string bgUnits = bkgWS->getAxis(0)->unit()->unitID();
  if (bgUnits != "TOF")
    throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() +
                                " should be in the units of TOF");

  if (!(bkgWS->getNumberHistograms() == 1 ||
        sourceWS->getNumberHistograms() == bkgWS->getNumberHistograms()))
    throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() +
                                " should have the same number of spectra as "
                                "source workspace or be a single histogram "
                                "workspace");

  auto WSUnit = sourceWS->getAxis(0)->unit();
  if (!WSUnit)
    throw std::invalid_argument(" Source Workspace: " + sourceWS->getName() +
                                " should have units");

  Geometry::IComponent_const_sptr source =
      sourceWS->getInstrument()->getSource();
  m_Sample = sourceWS->getInstrument()->getSample();
  if ((!source) || (!m_Sample))
    throw std::invalid_argument(
        "Instrument on Source workspace:" + sourceWS->getName() +
        "is not sufficiently defined: failed to get source and/or sample");
  m_L1 = source->getDistance(*m_Sample);

  // just in case.
  this->deleteUnitsConverters();
  // allocate the array of units converters to avoid units reallocation within a
  // loop
  m_WSUnit.assign(nThreads, NULL);
  for (int i = 0; i < nThreads; i++) {
    m_WSUnit[i] = WSUnit->clone();
  }

  m_singleValueBackground = false;
  if (bkgWS->getNumberHistograms() == 0)
    m_singleValueBackground = true;
  const MantidVec &dataX = bkgWS->dataX(0);
  const MantidVec &dataY = bkgWS->dataY(0);
  // const MantidVec& dataE = bkgWS->dataE(0);
  m_NBg = dataY[0];
  m_dtBg = dataX[1] - dataX[0];
  // m_ErrSq  = dataE[0]*dataE[0]; // needs further clarification

  m_Efix = this->getEi(sourceWS);
}

/** Performs a simple check to see if the sizes of two workspaces are compatible
 *for a binary operation
 *  In order to be size compatible then the larger workspace
 *  must divide be the size of the smaller workspace leaving no remainder
 *
 *  @param lhs :: the first workspace to compare
 *  @param rhs :: the second workspace to compare
 *  @retval "" The two workspaces are size compatible
 *  @retval "<reason why not compatible>" The two workspaces are NOT size
 *compatible
 */
std::string Divide::checkSizeCompatibility(
    const API::MatrixWorkspace_const_sptr lhs,
    const API::MatrixWorkspace_const_sptr rhs) const {
  // --- Check for event workspaces - different than workspaces 2D! ---

  // A SingleValueWorkspace on the right matches anything
  if (rhs->size() == 1)
    return "";

  // A SingleValueWorkspace on the left only matches if rhs was single value
  // too. Why are you using mantid to do simple math?!?
  if (lhs->size() == 1)
    return "The left side cannot contain a single value if the right side "
           "isn't also a single value.";

  // If RHS only has one value (1D vertical), the number of histograms needs to
  // match.
  // Each lhs spectrum will be divided by that scalar
  // std::cout << "rhs->blocksize() " << rhs->blocksize() << std::endl;
  // Are we allowing the division by different # of spectra, using detector IDs
  // to match up?
  if (m_AllowDifferentNumberSpectra ||
      (rhs->blocksize() == 1 &&
       lhs->getNumberHistograms() == rhs->getNumberHistograms())) {
    return "";
  }

  if (m_matchXSize) {
    // Past this point, for a 2D WS operation, we require the X arrays to match.
    // Note this only checks the first spectrum
    if (!WorkspaceHelpers::matchingBins(lhs, rhs, true)) {
      return "X arrays must match when dividing 2D workspaces.";
    }
  }

  // We don't need to check for matching bins for events. Yay events!
  const size_t rhsSpec = rhs->getNumberHistograms();

  // If the rhs has a single spectrum, then we can divide. The block size does
  // NOT need to match,
  if (rhsSpec == 1)
    return "";

  // Otherwise, the number of histograms needs to match, but the block size of
  // each does NOT need to match.

  if (lhs->getNumberHistograms() == rhs->getNumberHistograms()) {
    return "";
  } else {
    return "Number of histograms not identical.";
  }
}

/** Checks if workspaces input to Q1D or Qxy are reasonable
  @param dataWS data workspace
  @param binWS (WavelengthAdj) workpace that will be checked to see if it has one spectrum and the same number of bins as dataWS
  @param detectWS (PixelAdj) passing NULL for this wont raise an error, if set it will be checked this workspace has as many histograms as dataWS each with one bin
  @throw invalid_argument if the workspaces are not mututially compatible
*/
void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS, 
     API::MatrixWorkspace_const_sptr binAdj, API::MatrixWorkspace_const_sptr detectAdj)
{
  if ( dataWS->getNumberHistograms() < 1 )
  {
    throw std::invalid_argument("Empty data workspace passed, can not continue");
  }

  //it is not an error for these workspaces not to exist
  if (binAdj)
  {
    if ( binAdj->getNumberHistograms() != 1 )
    {
      throw std::invalid_argument("The WavelengthAdj workspace must have one spectrum");
    }
    if ( binAdj->readY(0).size() != dataWS->readY(0).size() )
    {
      throw std::invalid_argument("The WavelengthAdj workspace's bins must match those of the detector bank workspace");
    }
    MantidVec::const_iterator reqX = dataWS->readX(0).begin();
    MantidVec::const_iterator testX = binAdj->readX(0).begin();
    for ( ; reqX != dataWS->readX(0).end(); ++reqX, ++testX)
    {
      if ( *reqX != *testX )
      {
        throw std::invalid_argument("The WavelengthAdj workspace must have matching bins with the detector bank workspace");
      }
    }
    if ( binAdj->isDistribution() != dataWS->isDistribution() )
    {
      throw std::invalid_argument("The distrbution/raw counts status of the wavelengthAdj and DetBankWorkspace must be the same, use ConvertToDistribution");
    }
  }
  else if( ! dataWS->isDistribution() )
  {
    //throw std::invalid_argument("The data workspace must be a distrbution if there is no Wavelength dependent adjustment");
  }
  
  if (detectAdj)
  {
    if ( detectAdj->blocksize() != 1 )
    {
      throw std::invalid_argument("The PixelAdj workspace must point to a workspace with single bin spectra, as only the first bin is used");
    }
    if ( detectAdj->getNumberHistograms() != dataWS->getNumberHistograms() )
    {
      throw std::invalid_argument("The PixelAdj workspace must have one spectrum for each spectrum in the detector bank workspace");
    }
  }
}

/** Checks that the two input workspaces have non-overlapping spectra numbers and contributing detectors
 *  @param ws1 :: The first input workspace
 *  @param ws2 :: The second input workspace
 *  @param checkSpectra :: set to true to check for overlapping spectra numbers (non-sensical for event workspaces)
 *  @throw std::invalid_argument If there is some overlap
 */
void ConjoinWorkspaces::checkForOverlap(API::MatrixWorkspace_const_sptr ws1, API::MatrixWorkspace_const_sptr ws2, bool checkSpectra) const
{
  // make sure we should bother checking
  if (!this->getProperty("CheckOverlapping"))
    return;
  // Loop through the first workspace adding all the spectrum numbers & UDETS to a set
  std::set<specid_t> spectra;
  std::set<detid_t> detectors;
  const size_t& nhist1 = ws1->getNumberHistograms();
  for (size_t i = 0; i < nhist1; ++i)
  {
    const ISpectrum * spec = ws1->getSpectrum(i);
    const specid_t spectrum = spec->getSpectrumNo();
    spectra.insert(spectrum);
    const std::set<detid_t> & dets = spec->getDetectorIDs();
    std::set<detid_t>::const_iterator it;
    for (it = dets.begin(); it != dets.end(); ++it)
    {
      detectors.insert(*it);
    }
  }

  // Now go throught the spectrum numbers & UDETS in the 2nd workspace, making sure that there's no overlap
  const size_t& nhist2 = ws2->getNumberHistograms();
  for (size_t j = 0; j < nhist2; ++j)
  {
    const ISpectrum * spec = ws2->getSpectrum(j);
    const specid_t spectrum = spec->getSpectrumNo();
    if (checkSpectra)
    {
      if ( spectrum > 0 && spectra.find(spectrum) != spectra.end() )
      {
        g_log.error("The input workspaces have overlapping spectrum numbers");
        throw std::invalid_argument("The input workspaces have overlapping spectrum numbers");
      }
    }
    const std::set<detid_t> & dets = spec->getDetectorIDs();
    std::set<detid_t>::const_iterator it;
    for (it = dets.begin(); it != dets.end(); ++it)
    {
      if ( detectors.find(*it) != detectors.end() )
      {
        g_log.error("The input workspaces have common detectors");
        throw std::invalid_argument("The input workspaces have common detectors");
      }
    }
  }
}

/**
 * Creates the output workspace for this algorithm
 * @param inputWorkspace A parent workspace to initialize from.
 * @return A pointer to the output workspace.
 */
API::MatrixWorkspace_sptr Transpose::createOutputWorkspace(
    API::MatrixWorkspace_const_sptr inputWorkspace) {
  Mantid::API::Axis *yAxis = getVerticalAxis(inputWorkspace);
  const size_t oldNhist = inputWorkspace->getNumberHistograms();
  const auto &inX = inputWorkspace->x(0);
  const size_t oldYlength = inputWorkspace->blocksize();
  const size_t oldVerticalAxislength = yAxis->length();

  // The input Y axis may be binned so the new X data should be too
  size_t newNhist(oldYlength), newXsize(oldVerticalAxislength),
      newYsize(oldNhist);
  MatrixWorkspace_sptr outputWorkspace = inputWorkspace->cloneEmpty();
  outputWorkspace->initialize(newNhist, newXsize, newYsize);
  outputWorkspace->setTitle(inputWorkspace->getTitle());
  outputWorkspace->setComment(inputWorkspace->getComment());
  outputWorkspace->copyExperimentInfoFrom(inputWorkspace.get());
  outputWorkspace->setYUnit(inputWorkspace->YUnit());
  outputWorkspace->setYUnitLabel(inputWorkspace->YUnitLabel());
  outputWorkspace->setDistribution(inputWorkspace->isDistribution());

  // Create a new numeric axis for Y the same length as the old X array
  // Values come from input X
  API::NumericAxis *newYAxis(nullptr);
  if (inputWorkspace->isHistogramData()) {
    newYAxis = new API::BinEdgeAxis(inX.rawData());
  } else {
    newYAxis = new API::NumericAxis(inX.rawData());
  }

  newYAxis->unit() = inputWorkspace->getAxis(0)->unit();
  outputWorkspace->getAxis(0)->unit() = inputWorkspace->getAxis(1)->unit();
  outputWorkspace->replaceAxis(1, newYAxis);
  setProperty("OutputWorkspace", outputWorkspace);
  return outputWorkspace;
}

/** Calculates rebin parameters: the min and max bin boundaries and the
   logarithmic step. The aim is to have approx.
    the same number of bins as in the input workspace.
    @param workspace :: The workspace being rebinned
    @param min ::       (return) The calculated frame starting point
    @param max ::       (return) The calculated frame ending point
    @param step ::      (return) The calculated bin width
 */
void DiffractionFocussing::calculateRebinParams(
    const API::MatrixWorkspace_const_sptr &workspace, double &min, double &max,
    double &step) {

  min = std::numeric_limits<double>::max();
  // for min and max we need to iterate over the data block and investigate each
  // one
  int64_t length = workspace->getNumberHistograms();
  for (int64_t i = 0; i < length; i++) {
    auto &xVec = workspace->x(i);
    const double &localMin = xVec.front();
    const double &localMax = xVec.back();
    if (std::isfinite(localMin) && std::isfinite(localMax)) {
      min = std::min(min, localMin);
      max = std::max(max, localMax);
    }
  }

  if (min <= 0.)
    min = 1e-6;

  // step is easy
  double n = static_cast<double>(workspace->blocksize());
  step = (log(max) - log(min)) / n;
}

/** Method updates the column, which describes if current detector/spectra is
   masked
    It is used if one tries to process multiple workspaces obtained from a
   series of experiments  where the masked detectors can change */
void PreprocessDetectorsToMD::updateMasksState(
    const API::MatrixWorkspace_const_sptr &inputWS,
    DataObjects::TableWorkspace_sptr &targWS) {
  int *pMasksArray = targWS->getColDataArray<int>("detMask");
  if (!pMasksArray)
    throw std::invalid_argument(
        "target workspace " + targWS->getName() +
        " does not have defined masks column to update");

  size_t nHist = targWS->rowCount();
  const size_t nRows = inputWS->getNumberHistograms();
  if (nHist != nRows)
    throw std::invalid_argument(
        " source workspace " + inputWS->getName() + " and target workspace " +
        targWS->getName() +
        " are inconsistent as have different numner of detectors");

  uint32_t liveDetectorsCount(0);
  const auto &spectrumInfo = inputWS->spectrumInfo();
  for (size_t i = 0; i < nHist; i++) {
    if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
      continue;

    // if masked detectors state is not used, masked detectors just ignored;
    bool maskDetector = spectrumInfo.isMasked(i);
    *(pMasksArray + liveDetectorsCount) = maskDetector ? 1 : 0;

    liveDetectorsCount++;
  }
}

/** Create an output workspace of the appropriate (histogram or event) type and
 * copy over the data
 *  @param inputWS The input workspace
 */
API::MatrixWorkspace_sptr ConvertUnitsUsingDetectorTable::setupOutputWorkspace(
    const API::MatrixWorkspace_const_sptr inputWS) {
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");

  // If input and output workspaces are NOT the same, create a new workspace for
  // the output
  if (outputWS != inputWS) {
    if (m_inputEvents) {
      // Need to create by name as WorkspaceFactory otherwise spits out
      // Workspace2D when EventWS passed in
      outputWS = WorkspaceFactory::Instance().create(
          "EventWorkspace", inputWS->getNumberHistograms(), 2, 1);
      // Copy geometry etc. over
      WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS,
                                                        false);
      // Need to copy over the data as well
      EventWorkspace_const_sptr inputEventWS =
          boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
      boost::dynamic_pointer_cast<EventWorkspace>(outputWS)
          ->copyDataFrom(*inputEventWS);
    } else {
      // Create the output workspace
      outputWS = WorkspaceFactory::Instance().create(inputWS);
      // Copy the data over
      this->fillOutputHist(inputWS, outputWS);
    }
  }

  // Set the final unit that our output workspace will have
  outputWS->getAxis(0)->unit() = m_outputUnit;

  return outputWS;
}

/***
 * This will add the maximum spectrum number from ws1 to the data coming from ws2.
 *
 * @param ws1 The first workspace supplied to the algorithm.
 * @param ws2 The second workspace supplied to the algorithm.
 * @param output The workspace that is going to be returned by the algorithm.
 */
void ConjoinWorkspaces::fixSpectrumNumbers(API::MatrixWorkspace_const_sptr ws1, API::MatrixWorkspace_const_sptr /*ws2*/,
                                      API::MatrixWorkspace_sptr output)
{
  // make sure we should bother. If you don't check for overlap, you don't need to
  if (this->getProperty("CheckOverlapping"))
    return;

  // is everything possibly ok?
  specid_t min;
  specid_t max;
  getMinMax(output, min, max);
  if (max - min >= static_cast<specid_t>(output->getNumberHistograms())) // nothing to do then
    return;

  // information for remapping the spectra numbers
  specid_t ws1min;
  specid_t ws1max;
  getMinMax(ws1, ws1min, ws1max);

  // change the axis by adding the maximum existing spectrum number to the current value
  for (size_t i = ws1->getNumberHistograms(); i < output->getNumberHistograms(); i++)
  {
    specid_t origid;
    origid = output->getSpectrum(i)->getSpectrumNo();
    output->getSpectrum(i)->setSpectrumNo(origid + ws1max);
  }
  // To be deprecated:
  output->generateSpectraMap();
}

/**
 * Creates the output workspace for this algorithm
 * @param inputWorkspace A parent workspace to initialize from.
 * @return A pointer to the output workspace.
 */
API::MatrixWorkspace_sptr Transpose::createOutputWorkspace(
    API::MatrixWorkspace_const_sptr inputWorkspace) {
  Mantid::API::Axis *yAxis = getVerticalAxis(inputWorkspace);
  const size_t oldNhist = inputWorkspace->getNumberHistograms();
  const MantidVec &inX = inputWorkspace->readX(0);
  const size_t oldYlength = inputWorkspace->blocksize();
  const size_t oldVerticalAxislength = yAxis->length();

  // The input Y axis may be binned so the new X data should be too
  size_t newNhist(oldYlength), newXsize(oldVerticalAxislength),
      newYsize(oldNhist);
  MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
      inputWorkspace, newNhist, newXsize, newYsize);

  // Create a new numeric axis for Y the same length as the old X array
  // Values come from input X
  API::NumericAxis *newYAxis(nullptr);
  if (inputWorkspace->isHistogramData()) {
    newYAxis = new API::BinEdgeAxis(inX);
  } else {
    newYAxis = new API::NumericAxis(inX);
  }

  newYAxis->unit() = inputWorkspace->getAxis(0)->unit();
  outputWorkspace->getAxis(0)->unit() = inputWorkspace->getAxis(1)->unit();
  outputWorkspace->replaceAxis(1, newYAxis);
  setProperty("OutputWorkspace", outputWorkspace);
  return outputWorkspace;
}

/** Calculates rebin parameters: the min and max bin boundaries and the
   logarithmic step. The aim is to have approx.
    the same number of bins as in the input workspace.
    @param workspace :: The workspace being rebinned
    @param min ::       (return) The calculated frame starting point
    @param max ::       (return) The calculated frame ending point
    @param step ::      (return) The calculated bin width
 */
void DiffractionFocussing::calculateRebinParams(
    const API::MatrixWorkspace_const_sptr &workspace, double &min, double &max,
    double &step) {

  min = std::numeric_limits<double>::max();
  // for min and max we need to iterate over the data block and investigate each
  // one
  int64_t length = workspace->getNumberHistograms();
  for (int64_t i = 0; i < length; i++) {
    const MantidVec &xVec = workspace->readX(i);
    const double &localMin = xVec[0];
    const double &localMax = xVec[xVec.size() - 1];
    if (localMin != std::numeric_limits<double>::infinity() &&
        localMax != std::numeric_limits<double>::infinity()) {
      if (localMin < min)
        min = localMin;
      if (localMax > max)
        max = localMax;
    }
  }

  if (min <= 0.)
    min = 1e-6;

  // step is easy
  double n = static_cast<double>(workspace->blocksize());
  step = (log(max) - log(min)) / n;
}

  /**
   * Write bin masking information
   * @param ws :: The workspace
   * @return true for OK, false for error
   */
  bool NexusFileIO::writeNexusBinMasking(API::MatrixWorkspace_const_sptr ws) const
  {
    std::vector< int > spectra;
    std::vector< std::size_t > bins;
    std::vector< double > weights;
    int spectra_count = 0;
    int offset = 0;
    for(std::size_t i=0;i<ws->getNumberHistograms(); ++i)
    {
      if (ws->hasMaskedBins(i))
      {
        const API::MatrixWorkspace::MaskList& mList = ws->maskedBins(i);
        spectra.push_back(spectra_count);
        spectra.push_back(offset);
        API::MatrixWorkspace::MaskList::const_iterator it = mList.begin();
        for(;it != mList.end(); ++it)
        {
          bins.push_back(it->first);
          weights.push_back(it->second);
        }
        ++spectra_count;
        offset += static_cast<int>(mList.size());
      }
    }

    if (spectra_count == 0) return false;

    NXstatus status;

    // save spectra offsets as a 2d array of ints
    int dimensions[2];
    dimensions[0]=spectra_count;
    dimensions[1]=2;
    status=NXmakedata(fileID, "masked_spectra", NX_INT32, 2, dimensions);
    if(status==NX_ERROR) return false;
    NXopendata(fileID, "masked_spectra");
    const std::string description = "spectra index,offset in masked_bins and mask_weights";
    NXputattr(fileID, "description",  reinterpret_cast<void*>(const_cast<char*>(description.c_str())), static_cast<int>(description.size()+1), NX_CHAR);
    NXputdata(fileID, (void*)&spectra[0]);
    NXclosedata(fileID);

    // save masked bin indices
    dimensions[0]=static_cast<int>(bins.size());
    status=NXmakedata(fileID, "masked_bins", NX_INT32, 1, dimensions);
    if(status==NX_ERROR) return false;
    NXopendata(fileID, "masked_bins");
    NXputdata(fileID, (void*)&bins[0]);
    NXclosedata(fileID);

    // save masked bin weights
    dimensions[0]=static_cast<int>(bins.size());
    status=NXmakedata(fileID, "mask_weights", NX_FLOAT64, 1, dimensions);
    if(status==NX_ERROR) return false;
    NXopendata(fileID, "mask_weights");
    NXputdata(fileID, (void*)&weights[0]);
    NXclosedata(fileID);

    return true;
  }