/** Converts detector IDs to spectra indexes
*  @param WS :: the workspace on which the calculations are being performed
*  @param specNum1 :: spectrum number of the output of the first monitor
*  @param specNum2 :: spectrum number of the output of the second monitor
*  @return the indexes of the histograms created by the detector whose ID were passed
*  @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
*/
std::vector<size_t> GetEi::getMonitorSpecIndexs(API::MatrixWorkspace_const_sptr WS, specid_t specNum1, specid_t specNum2) const
{// getting spectra numbers from detector IDs is hard because the map works the other way, getting index numbers from spectra numbers has the same problem and we are about to do both
  std::vector<size_t> specInds;
  
  // get the index number of the histogram for the first monitor
  std::vector<specid_t> specNumTemp(&specNum1, &specNum1+1);
  WS->getIndicesFromSpectra(specNumTemp, specInds);
  if ( specInds.size() != 1 )
  {// the monitor spectrum isn't present in the workspace, we can't continue from here
    g_log.error() << "Couldn't find the first monitor spectrum, number " << specNum1 << std::endl;
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum1);
  }

  // nowe the second monitor
  std::vector<size_t> specIndexTemp;
  specNumTemp[0] = specNum2;
  WS->getIndicesFromSpectra(specNumTemp, specIndexTemp);
  if ( specIndexTemp.size() != 1 )
  {// the monitor spectrum isn't present in the workspace, we can't continue from here
    g_log.error() << "Couldn't find the second monitor spectrum, number " << specNum2 << std::endl;
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum2);
  }
  
  specInds.push_back(specIndexTemp[0]);
  return specInds;
}

/** 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;
}

/** Calculate the normalization term for each output bin
*  @param wavStart [in] the index number of the first bin in the input
* wavelengths that is actually being used
*  @param specInd [in] the spectrum to calculate
*  @param pixelAdj [in] if not NULL this is workspace contains single bins with
* the adjustments, e.g. detector efficencies, for the given spectrum index
*  @param wavePixelAdj [in] if not NULL this is workspace that contains the
* adjustments for the pixels and wavelenght dependend values.
*  @param binNorms [in] pointer to a contigious array of doubles that are the
* wavelength correction from waveAdj workspace, can be NULL
*  @param binNormEs [in] pointer to a contigious array of doubles which
* corrospond to the corrections and are their errors, can be NULL
*  @param norm [out] normalization for each bin, including soild angle, pixel
* correction, the proportion that is not masked and the normalization workspace
*  @param normETo2 [out] this pointer must point to the end of the norm array,
* it will be filled with the total of the error on the normalization
*/
void Q1D2::calculateNormalization(const size_t wavStart, const size_t specInd,
                                  API::MatrixWorkspace_const_sptr pixelAdj,
                                  API::MatrixWorkspace_const_sptr wavePixelAdj,
                                  double const *const binNorms,
                                  double const *const binNormEs,
                                  const MantidVec::iterator norm,
                                  const MantidVec::iterator normETo2) const {
  double detectorAdj, detAdjErr;
  pixelWeight(pixelAdj, specInd, detectorAdj, detAdjErr);
  // use that the normalization array ends at the start of the error array
  for (MantidVec::iterator n = norm, e = normETo2; n != normETo2; ++n, ++e) {
    *n = detectorAdj;
    *e = detAdjErr *detAdjErr;
  }

  if (binNorms && binNormEs) {
    if (wavePixelAdj)
      // pass the iterator for the wave pixel Adj dependent
      addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2,
                 wavePixelAdj->readY(specInd).begin() + wavStart,
                 wavePixelAdj->readE(specInd).begin() + wavStart);
    else
      addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2);
  }
  normToMask(wavStart, specInd, norm, normETo2);
}

/** Get the time vector for the TimeSeriesProperty to which the entries is to
 * set
 * @brief AddSampleLog::getTimes
 * @param dataws
 * @param workspace_index
 * @param is_epoch
 * @param is_second
 * @param run_obj
 * @return
 */
std::vector<Types::Core::DateAndTime>
AddSampleLog::getTimes(API::MatrixWorkspace_const_sptr dataws,
                       int workspace_index, bool is_epoch, bool is_second,
                       API::Run &run_obj) {
  // get run start time
  int64_t timeshift(0);
  if (!is_epoch) {
    // get the run start time
    Types::Core::DateAndTime run_start_time = getRunStart(run_obj);
    timeshift = run_start_time.totalNanoseconds();
  }

  // set up the time vector
  std::vector<Types::Core::DateAndTime> timevec;
  size_t vecsize = dataws->readX(workspace_index).size();
  for (size_t i = 0; i < vecsize; ++i) {
    double timedbl = dataws->readX(workspace_index)[i];
    if (is_second)
      timedbl *= 1.E9;
    int64_t entry_i64 = static_cast<int64_t>(timedbl);
    Types::Core::DateAndTime entry(timeshift + entry_i64);
    timevec.push_back(entry);
  }

  return timevec;
}

 /**
  * Checks if the spectra at the given index of either input workspace is masked. If so then the output spectra has zeroed data
  * and is also masked. 
  * @param lhs :: A pointer to the left-hand operand
  * @param rhs :: A pointer to the right-hand operand
  * @param index :: The workspace index to check
  * @param out :: A pointer to the output workspace
  * @returns True if further processing is not required on the spectra, false if the binary operation should be performed.
  */
 bool BinaryOperation::propagateSpectraMask(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs, 
     const int64_t index, API::MatrixWorkspace_sptr out)
 {
   bool continueOp(true);
   IDetector_const_sptr det_lhs, det_rhs;
   try
   {
     det_lhs = lhs->getDetector(index);
     det_rhs = rhs->getDetector(index);
   }
   catch(std::runtime_error &)
   {
   }
   catch(std::domain_error &)
   {
     // try statement will throw a domain_error when the axis is not a spectra axis.
     return continueOp;
   }
   if( (det_lhs && det_lhs->isMasked()) || ( det_rhs && det_rhs->isMasked()) )
   {
     continueOp = false;
     out->maskWorkspaceIndex(index);
   }
   return continueOp;
 }

/** 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++;
  }
}

/** Calculates the parameters to hand to the Rebin algorithm. Specifies the new
 * binning, bin-by-bin,
 *  to cover the full range covered by the two input workspaces. In regions of
 * overlap, the bins from
 *  the workspace having the wider bins are taken. Note that because the list of
 * input workspaces
 *  is sorted, ws1 will always start before (or at the same point as) ws2.
 *  @param ws1 ::    The first input workspace. Will start before ws2.
 *  @param ws2 ::    The second input workspace.
 *  @param params :: A reference to the vector of rebinning parameters
 */
void MergeRuns::calculateRebinParams(const API::MatrixWorkspace_const_sptr &ws1,
                                     const API::MatrixWorkspace_const_sptr &ws2,
                                     std::vector<double> &params) const {
    const MantidVec &X1 = ws1->readX(0);
    const MantidVec &X2 = ws2->readX(0);
    const double end1 = X1.back();
    const double start2 = X2.front();
    const double end2 = X2.back();

    if (end1 <= start2) {
        // First case is if there's no overlap between the workspaces
        this->noOverlapParams(X1, X2, params);
    } else {
        // Add the bins from the first workspace up to the start of the overlap
        params.push_back(X1[0]);
        int64_t i;
        for (i = 1; X1[i] <= start2; ++i) {
            params.push_back(X1[i] - X1[i - 1]);
            params.push_back(X1[i]);
        }
        // If the range of workspace2 is completely within that of workspace1, then
        // call the
        // 'inclusion' routine. Otherwise call the standard 'intersection' one.
        if (end1 < end2) {
            this->intersectionParams(X1, i, X2, params);
        } else {
            this->inclusionParams(X1, i, X2, params);
        }
    }
}

/** Method returns the efixed or Ei value stored in properties of the input
 *workspace.
 *  Indirect instruments can have eFxed and Direct instruments can have Ei
 *defined as the properties of the workspace.
 *
 *  This method provide guess for efixed for all other kind of instruments.
 *Correct indirect instrument will overwrite
 *  this value while wrongly defined or different types of instruments will
 *provide the value of "Ei" property (log value)
 *  or undefined if "Ei" property is not found.
 *
 */
double PreprocessDetectorsToMD::getEi(
    const API::MatrixWorkspace_const_sptr &inputWS) const {
  double Efi = std::numeric_limits<double>::quiet_NaN();

  // is Ei on workspace properties? (it can be defined for some reason if
  // detectors do not have one, and then it would exist as Ei)
  bool EiFound(false);
  try {
    Efi = inputWS->run().getPropertyValueAsType<double>("Ei");
    EiFound = true;
  } catch (Kernel::Exception::NotFoundError &) {
  }
  // try to get Efixed as property on a workspace, obtained for indirect
  // instrument
  bool eFixedFound(false);
  if (!EiFound) {
    try {
      Efi = inputWS->run().getPropertyValueAsType<double>("eFixed");
      eFixedFound = true;
    } catch (Kernel::Exception::NotFoundError &) {
    }
  }

  if (!(EiFound || eFixedFound))
    g_log.debug() << " Ei/eFixed requested but have not been found\n";

  return Efi;
}

    /**
     * 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;
    }

/** 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;
}

size_t
EnggDiffFittingModel::guessBankID(API::MatrixWorkspace_const_sptr ws) const {
  const static std::string bankIDName = "bankid";
  if (ws->run().hasProperty(bankIDName)) {
    const auto log = dynamic_cast<Kernel::PropertyWithValue<int> *>(
        ws->run().getLogData(bankIDName));
    return boost::lexical_cast<size_t>(log->value());
  }

  // couldn't get it from sample logs - try using the old naming convention
  const std::string name = ws->getName();
  std::vector<std::string> chunks;
  boost::split(chunks, name, boost::is_any_of("_"));
  bool isNum = isDigit(chunks.back());
  if (!chunks.empty() && isNum) {
    try {
      return boost::lexical_cast<size_t>(chunks.back());
    } catch (boost::exception &) {
      // If we get a bad cast or something goes wrong then
      // the file is probably not what we were expecting
      // so throw a runtime error
      throw std::runtime_error(
          "Failed to fit file: The data was not what is expected. "
          "Does the file contain a focused workspace?");
    }
  }

  throw std::runtime_error("Could not guess run number from input workspace. "
                           "Are you sure it has been focused correctly?");
}

/** 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.");
      }
    }
  }
}

/** Execute the algorithm.
 */
void DampSq::exec()
{
    // TODO Auto-generated execute stub

    // 1. Generate new workspace
    API::MatrixWorkspace_const_sptr isqspace = getProperty("InputWorkspace");
    API::MatrixWorkspace_sptr osqspace = WorkspaceFactory::Instance().create(isqspace, 1, isqspace->size(), isqspace->size());

    int mode = getProperty("Mode");
    double qmax = getProperty("QMax");

    if (mode < 1 || mode > 4) {
        g_log.error("Damp mode can only be 1, 2, 3, or 4");
        return;
    }

    // 2. Get access to all
    const MantidVec& iQVec = isqspace->dataX(0);
    const MantidVec& iSVec = isqspace->dataY(0);
    const MantidVec& iEVec = isqspace->dataE(0);

    MantidVec& oQVec = osqspace->dataX(0);
    MantidVec& oSVec = osqspace->dataY(0);
    MantidVec& oEVec = osqspace->dataE(0);

    // 3. Calculation
    double dqmax = qmax - iQVec[0];

    double damp;
    for (unsigned int i = 0; i < iQVec.size(); i ++) {
        // a) calculate damp coefficient
        switch (mode) {
        case 1:
            damp = dampcoeff1(iQVec[i], qmax, dqmax);
            break;
        case 2:
            damp = dampcoeff2(iQVec[i], qmax, dqmax);;
            break;
        case 3:
            damp = dampcoeff3(iQVec[i], qmax, dqmax);;
            break;
        case 4:
            damp = dampcoeff4(iQVec[i], qmax, dqmax);;
            break;
        default:
            damp = 0;
            break;
        }
        // b) calculate new S(q)
        oQVec[i] = iQVec[i];
        oSVec[i] = 1 + damp*(iSVec[i]-1);
        oEVec[i] = damp*iEVec[i];
    }  // i

    // 4. Over
    setProperty("OutputWorkspace", osqspace);

    return;
}

  /**
   * 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;
  }

// read the monitors list from the workspace and try to do it once for any
// particular ws;
bool MonIDPropChanger::monitorIdReader(
    API::MatrixWorkspace_const_sptr inputWS) const {
  // no workspace
  if (!inputWS)
    return false;

  // no instrument
  Geometry::Instrument_const_sptr pInstr = inputWS->getInstrument();
  if (!pInstr)
    return false;

  std::vector<detid_t> mon = pInstr->getMonitors();
  if (mon.empty()) {
    if (iExistingAllowedValues.empty()) {
      return false;
    } else {
      iExistingAllowedValues.clear();
      return true;
    }
  }
  // are these monitors really there?
  // got the index of correspondent spectra.
  std::vector<size_t> indexList = inputWS->getIndicesFromDetectorIDs(mon);
  if (indexList.empty()) {
    if (iExistingAllowedValues.empty()) {
      return false;
    } else {
      iExistingAllowedValues.clear();
      return true;
    }
  }
  // index list can be less or equal to the mon list size (some monitors do not
  // have spectra)
  size_t mon_count =
      (mon.size() < indexList.size()) ? mon.size() : indexList.size();
  std::vector<int> allowed_values(mon_count);
  for (size_t i = 0; i < mon_count; i++) {
    allowed_values[i] = mon[i];
  }

  // are known values the same as the values we have just identified?
  if (iExistingAllowedValues.size() != mon_count) {
    iExistingAllowedValues.clear();
    iExistingAllowedValues.assign(allowed_values.begin(), allowed_values.end());
    return true;
  }
  // the monitor list has the same size as before. Is it equivalent to the
  // existing one?
  bool values_redefined = false;
  for (size_t i = 0; i < mon_count; i++) {
    if (iExistingAllowedValues[i] != allowed_values[i]) {
      values_redefined = true;
      iExistingAllowedValues[i] = allowed_values[i];
    }
  }
  return values_redefined;
}

/** Get the values (in integer) of the TimeSeriesProperty's entries from input
 * workspace
 * @brief AddSampleLog::getIntValues
 * @param dataws
 * @param workspace_index
 * @return
 */
std::vector<int>
AddSampleLog::getIntValues(API::MatrixWorkspace_const_sptr dataws,
                           int workspace_index) {
  std::vector<int> valuevec;
  size_t vecsize = dataws->readY(workspace_index).size();
  for (size_t i = 0; i < vecsize; ++i)
    valuevec.push_back(static_cast<int>(dataws->readY(workspace_index)[i]));

  return valuevec;
}

    /**
     * Create a masking workspace to return.
     *
     * @param inputWS The workspace to initialize from. The instrument is copied from this.
     */
    DataObjects::MaskWorkspace_sptr DetectorDiagnostic::generateEmptyMask(API::MatrixWorkspace_const_sptr inputWS)
    {
      // Create a new workspace for the results, copy from the input to ensure that we copy over the instrument and current masking
      DataObjects::MaskWorkspace_sptr maskWS(new DataObjects::MaskWorkspace());
      maskWS->initialize(inputWS->getNumberHistograms(), 1, 1);
      WorkspaceFactory::Instance().initializeFromParent(inputWS, maskWS, false);
      maskWS->setTitle(inputWS->getTitle());

      return maskWS;
    }

    /**
     * A map detector ID and Q ranges
     * This method looks unnecessary as it could be calculated on the fly but
     * the parallelization means that lazy instantation slows it down due to the
     * necessary CRITICAL sections required to update the cache. The Q range
     * values are required very frequently so the total time is more than
     * offset by this precaching step
     */
    void SofQW2::initThetaCache(API::MatrixWorkspace_const_sptr workspace)
    {
      const size_t nhist = workspace->getNumberHistograms();
      m_thetaPts = std::vector<double>(nhist);
      size_t ndets(0);
      double minTheta(DBL_MAX), maxTheta(-DBL_MAX);

      for(int64_t i = 0 ; i < (int64_t)nhist; ++i) //signed for OpenMP
      {

        m_progress->report("Calculating detector angles");
        IDetector_const_sptr det;
        try
        {
          det = workspace->getDetector(i);
          // Check to see if there is an EFixed, if not skip it
          try
          {
            m_EmodeProperties.getEFixed(det);
          }
          catch(std::runtime_error&)
          {
            det.reset();
          }
        }
        catch(Kernel::Exception::NotFoundError&)
        {
          // Catch if no detector. Next line tests whether this happened - test placed
          // outside here because Mac Intel compiler doesn't like 'continue' in a catch
          // in an openmp block.
        }
        // If no detector found, skip onto the next spectrum
        if( !det || det->isMonitor() )
        {
          m_thetaPts[i] = -1.0; // Indicates a detector to skip
        }
        else
        {
          ++ndets;
          const double theta = workspace->detectorTwoTheta(det);
          m_thetaPts[i] = theta;
          if( theta < minTheta )
          {
            minTheta = theta;
          }
          else if( theta > maxTheta )
          {
            maxTheta = theta;
          }
        }
      }

      m_thetaWidth = (maxTheta - minTheta)/static_cast<double>(ndets);
      g_log.information() << "Calculated detector width in theta=" << (m_thetaWidth*180.0/M_PI) << " degrees.\n";
    }

/** Calculate the derivatives of the given order from the interpolated points
 *
 * @param inputWorkspace :: The input workspace
 * @param outputWorkspace :: The output workspace
 * @param order :: The order of derivatives to calculate
 */
void SplineInterpolation::calculateDerivatives(
    API::MatrixWorkspace_const_sptr inputWorkspace,
    API::MatrixWorkspace_sptr outputWorkspace, int order) const {
  // get x and y parameters from workspaces
  const size_t nData = inputWorkspace->y(0).size();
  const double *xValues = &(inputWorkspace->x(0)[0]);
  double *yValues = &(outputWorkspace->mutableY(order - 1)[0]);

  // calculate the derivatives
  m_cspline->derivative1D(yValues, xValues, nData, order);
}

T ConvertEmptyToTof::getPropertyFromRun(API::MatrixWorkspace_const_sptr inputWS,
                                        const std::string &propertyName) {
  if (inputWS->run().hasProperty(propertyName)) {
    Kernel::Property *prop = inputWS->run().getProperty(propertyName);
    return boost::lexical_cast<T>(prop->value());
  } else {
    std::string mesg =
        "No '" + propertyName + "' property found in the input workspace....";
    throw std::runtime_error(mesg);
  }
}