//----------------------------------------------------------------------------------------------
void LoadPDFgetNFile::setUnit(Workspace2D_sptr ws) {
  // 1. Set X
  string xcolname = mColumnNames[0];

  if (xcolname.compare("Q") == 0) {
    string unit = "MomentumTransfer";
    ws->getAxis(0)->setUnit(unit);
  } else if (xcolname.compare("r") == 0) {
    ws->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
    Unit_sptr unit = ws->getAxis(0)->unit();
    boost::shared_ptr<Units::Label> label =
        boost::dynamic_pointer_cast<Units::Label>(unit);
    label->setLabel("AtomicDistance", "Angstrom");
  } else {
    stringstream errss;
    errss << "X axis " << xcolname << " is not supported for unit. " << endl;
    g_log.warning() << errss.str() << endl;
  }

  // 2. Set Y
  string ycolname = mColumnNames[1];
  string ylabel("");
  if (ycolname.compare("G(r)") == 0) {
    ylabel = "PDF";
  } else if (ycolname.compare("S") == 0) {
    ylabel = "S";
  } else {
    ylabel = "Intensity";
  }
  ws->setYUnitLabel(ylabel);

  return;
}

void MDHistoToWorkspace2D::recurseData(IMDHistoWorkspace_sptr inWS,
                                       Workspace2D_sptr outWS,
                                       size_t currentDim, coord_t *pos) {
  boost::shared_ptr<const IMDDimension> dim = inWS->getDimension(currentDim);
  if (currentDim == rank - 1) {
    MantidVec &Y = outWS->dataY(currentSpectra);
    for (unsigned int j = 0; j < dim->getNBins(); j++) {
      pos[currentDim] = dim->getX(j);
      Y[j] = inWS->getSignalAtCoord(
          pos, static_cast<Mantid::API::MDNormalization>(0));
    }
    MantidVec &E = outWS->dataE(currentSpectra);
    // MSVC compiler can't figure out the correct overload with out the function
    // cast on sqrt
    std::transform(Y.begin(), Y.end(), E.begin(),
                   (double (*)(double))std::sqrt);
    std::vector<double> xData;
    for (unsigned int i = 0; i < dim->getNBins(); i++) {
      xData.push_back(dim->getX(i));
    }
    outWS->setX(currentSpectra, xData);
    outWS->getSpectrum(currentSpectra)
        ->setSpectrumNo(static_cast<specid_t>(currentSpectra));
    currentSpectra++;
  } else {
    // recurse deeper
    for (int i = 0; i < static_cast<int>(dim->getNBins()); i++) {
      pos[currentDim] = dim->getX(i);
      recurseData(inWS, outWS, currentDim + 1, pos);
    }
  }
}

/** Remove peaks from a input workspace
  */
Workspace2D_sptr
RemovePeaks::removePeaks(API::MatrixWorkspace_const_sptr dataws, int wsindex,
                         double numfwhm) {
  // Check
  if (m_vecPeakCentre.empty())
    throw runtime_error("RemovePeaks has not been setup yet. ");

  // Initialize vectors
  const MantidVec &vecX = dataws->readX(wsindex);
  const MantidVec &vecY = dataws->readY(wsindex);
  const MantidVec &vecE = dataws->readE(wsindex);

  size_t sizex = vecX.size();
  vector<bool> vec_useX(sizex, true);

  // Exclude regions
  size_t numbkgdpoints =
      excludePeaks(vecX, vec_useX, m_vecPeakCentre, m_vecPeakFWHM, numfwhm);
  size_t numbkgdpointsy = numbkgdpoints;
  size_t sizey = vecY.size();
  if (sizex > sizey)
    --numbkgdpointsy;

  // Construct output workspace
  Workspace2D_sptr outws = boost::dynamic_pointer_cast<Workspace2D>(
      WorkspaceFactory::Instance().create("Workspace2D", 1, numbkgdpoints,
                                          numbkgdpointsy));
  outws->getAxis(0)->setUnit(dataws->getAxis(0)->unit()->unitID());
  MantidVec &outX = outws->dataX(0);
  MantidVec &outY = outws->dataY(0);
  MantidVec &outE = outws->dataE(0);
  size_t index = 0;
  for (size_t i = 0; i < sizex; ++i) {
    if (vec_useX[i]) {
      if (index >= numbkgdpoints)
        throw runtime_error("Programming logic error (1)");
      outX[index] = vecX[i];
      ++index;
    }
  }
  index = 0;
  for (size_t i = 0; i < sizey; ++i) {
    if (vec_useX[i]) {
      if (index >= numbkgdpointsy)
        throw runtime_error("Programming logic error (2)");
      outY[index] = vecY[i];
      outE[index] = vecE[i];
      ++index;
    }
  }

  return outws;
}

  /** Fit function
    * Minimizer: "Levenberg-MarquardtMD"/"Simplex"
   */
  bool RefinePowderInstrumentParameters2::doFitFunction(IFunction_sptr function, Workspace2D_sptr dataws, int wsindex,
                                                        string minimizer, int numiters, double& chi2, string& fitstatus)
  {
    // 0. Debug output
    stringstream outss;
    outss << "Fit function: " << m_positionFunc->asString() << endl << "Data To Fit: \n";
    for (size_t i = 0; i < dataws->readX(0).size(); ++i)
      outss << dataws->readX(wsindex)[i] << "\t\t" << dataws->readY(wsindex)[i] << "\t\t"
            << dataws->readE(wsindex)[i] << "\n";
    g_log.information() << outss.str();

    // 1. Create and setup fit algorithm
    API::IAlgorithm_sptr fitalg = createChildAlgorithm("Fit", 0.0, 0.2, true);
    fitalg->initialize();

    fitalg->setProperty("Function", function);
    fitalg->setProperty("InputWorkspace", dataws);
    fitalg->setProperty("WorkspaceIndex", wsindex);
    fitalg->setProperty("Minimizer", minimizer);
    fitalg->setProperty("CostFunction", "Least squares");
    fitalg->setProperty("MaxIterations", numiters);
    fitalg->setProperty("CalcErrors", true);

    // 2. Fit
    bool successfulfit = fitalg->execute();
    if (!fitalg->isExecuted() || ! successfulfit)
    {
      // Early return due to bad fit
      g_log.warning("Fitting to instrument geometry function failed. ");
      chi2 = DBL_MAX;
      fitstatus = "Minimizer throws exception.";
      return false;
    }

    // 3. Understand solution
    chi2 = fitalg->getProperty("OutputChi2overDoF");
    string tempfitstatus = fitalg->getProperty("OutputStatus");
    fitstatus = tempfitstatus;

    bool goodfit = fitstatus.compare("success") == 0;

    stringstream dbss;
    dbss << "Fit Result (GSL):  Chi^2 = " << chi2
         << "; Fit Status = " << fitstatus << ", Return Bool = " << goodfit << std::endl;
    vector<string> funcparnames = function->getParameterNames();
    for (size_t i = 0; i < funcparnames.size(); ++i)
      dbss << funcparnames[i] << " = " << setw(20) << function->getParameter(funcparnames[i])
           << " +/- " << function->getError(i) << "\n";
    g_log.debug() << dbss.str();

    return goodfit;
  }

/**
 * Reads the data (FITS matrix) from a single FITS file into a
 * workspace (directly into the spectra, using one spectrum per image
 * row).
 *
 * @param fileInfo information on the FITS file to load, including its path
 * @param cmpp centimeters per pixel, to scale/normalize values
 * @param ws workspace with the required dimensions
 * @param buffer pre-allocated buffer to read from file
 *
 * @throws std::runtime_error if there are file input issues
 */
void LoadFITS::readDataToWorkspace(const FITSInfo &fileInfo, double cmpp,
                                   Workspace2D_sptr ws,
                                   std::vector<char> &buffer) {
  const size_t bytespp = (fileInfo.bitsPerPixel / 8);
  const size_t len = m_pixelCount * bytespp;
  readInBuffer(fileInfo, buffer, len);

  const size_t nrows(fileInfo.axisPixelLengths[1]),
      ncols(fileInfo.axisPixelLengths[0]);
  // Treat buffer as a series of bytes
  uint8_t *buffer8 = reinterpret_cast<uint8_t *>(buffer.data());

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int i = 0; i < static_cast<int>(nrows); ++i) {
    auto &xVals = ws->mutableX(i);
    auto &yVals = ws->mutableY(i);
    auto &eVals = ws->mutableE(i);
    xVals = static_cast<double>(i) * cmpp;

    for (size_t j = 0; j < ncols; ++j) {
      // Map from 2D->1D index
      const size_t start = ((i * (bytespp)) * nrows) + (j * (bytespp));
      uint8_t const *const buffer8Start = buffer8 + start;
      // Reverse byte order of current value. Make sure we allocate enough
      // enough space to hold the size
      uint8_t byteValue[g_maxBytesPP];
      std::reverse_copy(buffer8Start, buffer8Start + bytespp, byteValue);

      double val = 0;
      if (fileInfo.bitsPerPixel == 8) {
        val = toDouble<uint8_t>(byteValue);
      } else if (fileInfo.bitsPerPixel == 16) {
        val = toDouble<uint16_t>(byteValue);
      } else if (fileInfo.bitsPerPixel == 32 && !fileInfo.isFloat) {
        val = toDouble<uint32_t>(byteValue);
      } else if (fileInfo.bitsPerPixel == 64 && !fileInfo.isFloat) {
        val = toDouble<uint32_t>(byteValue);
      } else if (fileInfo.bitsPerPixel == 32 && fileInfo.isFloat) {
        val = toDouble<float>(byteValue);
      } else if (fileInfo.bitsPerPixel == 64 && fileInfo.isFloat) {
        val = toDouble<double>(byteValue);
      }

      val = fileInfo.scale * val - fileInfo.offset;
      yVals[j] = val;
      eVals[j] = sqrt(val);
    }
  }
}

/** Sets error of workspace to specified value
  *
  * Since an estimation of the error is calculated from background counts, this
  *value is assigned to the workspace via this method.
  *
  * @param correlationWorkspace :: Workspace containing the correlation spectrum
  *on which the peak search was performed.
  * @param error :: Error that is set on the workspace.
  */
void
PoldiPeakSearch::setErrorsOnWorkspace(Workspace2D_sptr correlationWorkspace,
                                      double error) const {
  MantidVec &errors = correlationWorkspace->dataE(0);

  std::fill(errors.begin(), errors.end(), error);
}

/*
 * Generate a SANS test workspace, with instrument geometry.
 * The geometry is the SANSTEST geometry, with a 30x30 pixel 2D detector.
 *
 * @param workspace: name of the workspace to be created.
 */
Workspace2D_sptr SANSInstrumentCreationHelper::createSANSInstrumentWorkspace(
    std::string workspace) {
  // Create a test workspace with test data with a well defined peak
  // The test instrument has two monitor channels
  Workspace2D_sptr ws = WorkspaceCreationHelper::create2DWorkspace123(
      nBins * nBins + nMonitors, 1, 1);
  AnalysisDataService::Instance().addOrReplace(workspace, ws);
  ws->getAxis(0)->unit() =
      Mantid::Kernel::UnitFactory::Instance().create("Wavelength");
  ws->setYUnit("");

  // Load instrument geometry
  runLoadInstrument("SANSTEST", ws);
  runLoadMappingTable(ws, nBins, nBins);

  return ws;
}

/**
 * Writes a single workspace into the file
 * @param workspace the workspace to get data from
 * @param nxFile the nexus file to save data into
 */
void SaveNXTomo::writeSingleWorkspace(const Workspace2D_sptr workspace,
                                      ::NeXus::File &nxFile) {
  try {
    nxFile.openPath("/entry1/tomo_entry/data");
  } catch (...) {
    throw std::runtime_error("Unable to create a valid NXTomo file");
  }

  int numFiles = 0;
  nxFile.getAttr<int>("NumFiles", numFiles);

  // Change slab start to after last data position
  m_slabStart[0] = numFiles;
  m_slabSize[0] = 1;

  // Set the rotation value for this WS
  std::vector<double> rotValue;
  rotValue.push_back(0);

  if (workspace->run().hasProperty("Rotation")) {
    std::string tmpVal = workspace->run().getLogData("Rotation")->value();
    try {
      rotValue[0] = boost::lexical_cast<double>(tmpVal);
    } catch (...) {
    }
    // Invalid Cast is handled below
  }

  nxFile.openData("rotation_angle");
  nxFile.putSlab(rotValue, numFiles, 1);
  nxFile.closeData();

  // Copy data out, remake data with dimension of old size plus new elements.
  // Insert previous data.
  nxFile.openData("data");

  double *dataArr = new double[m_spectraCount];

  for (int64_t i = 0; i < m_dimensions[1]; ++i) {
    for (int64_t j = 0; j < m_dimensions[2]; ++j) {
      dataArr[i * m_dimensions[1] + j] =
          workspace->dataY(i * m_dimensions[1] + j)[0];
    }
  }

  nxFile.putSlab(dataArr, m_slabStart, m_slabSize);

  nxFile.closeData();

  nxFile.putAttr("NumFiles", numFiles + 1);

  nxFile.closeGroup();

  // Write additional log information, intensity and image key
  writeLogValues(workspace, nxFile, numFiles);
  writeIntensityValue(workspace, nxFile, numFiles);
  writeImageKeyValue(workspace, nxFile, numFiles);

  ++numFiles;

  delete[] dataArr;
}

/** Filter non-background data points out and create a background workspace
  */
Workspace2D_sptr
ProcessBackground::filterForBackground(BackgroundFunction_sptr bkgdfunction) {
  double posnoisetolerance = getProperty("NoiseTolerance");
  double negnoisetolerance = getProperty("NegativeNoiseTolerance");
  if (isEmpty(negnoisetolerance))
    negnoisetolerance = posnoisetolerance;

  // Calcualte theoretical values
  const std::vector<double> x = m_dataWS->readX(m_wsIndex);
  API::FunctionDomain1DVector domain(x);
  API::FunctionValues values(domain);
  bkgdfunction->function(domain, values);

  g_log.information() << "Function used to select background points : "
                      << bkgdfunction->asString() << "\n";

  // Optional output
  string userbkgdwsname = getPropertyValue("UserBackgroundWorkspace");
  if (userbkgdwsname.size() == 0)
    throw runtime_error("In mode SelectBackgroundPoints, "
                        "UserBackgroundWorkspace must be given!");

  size_t sizex = domain.size();
  size_t sizey = values.size();
  MatrixWorkspace_sptr visualws = boost::dynamic_pointer_cast<MatrixWorkspace>(
      WorkspaceFactory::Instance().create("Workspace2D", 4, sizex, sizey));
  for (size_t i = 0; i < sizex; ++i) {
    for (size_t j = 0; j < 4; ++j) {
      visualws->dataX(j)[i] = domain[i];
    }
  }
  for (size_t i = 0; i < sizey; ++i) {
    visualws->dataY(0)[i] = values[i];
    visualws->dataY(1)[i] = m_dataWS->readY(m_wsIndex)[i] - values[i];
    visualws->dataY(2)[i] = posnoisetolerance;
    visualws->dataY(3)[i] = -negnoisetolerance;
  }
  setProperty("UserBackgroundWorkspace", visualws);

  // Filter for background
  std::vector<double> vecx, vecy, vece;
  for (size_t i = 0; i < domain.size(); ++i) {
    // double y = m_dataWS->readY(m_wsIndex)[i];
    // double theoryy = values[i]; y-theoryy
    double purey = visualws->readY(1)[i];
    if (purey < posnoisetolerance && purey > -negnoisetolerance) {
      // Selected
      double x = domain[i];
      double y = m_dataWS->readY(m_wsIndex)[i];
      double e = m_dataWS->readE(m_wsIndex)[i];
      vecx.push_back(x);
      vecy.push_back(y);
      vece.push_back(e);
    }
  }
  g_log.information() << "Found " << vecx.size() << " background points out of "
                      << m_dataWS->readX(m_wsIndex).size()
                      << " total data points. "
                      << "\n";

  // Build new workspace for OutputWorkspace
  size_t nspec = 3;
  Workspace2D_sptr outws =
      boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", nspec,
                                                   vecx.size(), vecy.size()));
  for (size_t i = 0; i < vecx.size(); ++i) {
    for (size_t j = 0; j < nspec; ++j)
      outws->dataX(j)[i] = vecx[i];
    outws->dataY(0)[i] = vecy[i];
    outws->dataE(0)[i] = vece[i];
  }

  return outws;
}

/** Select background automatically
 */
DataObjects::Workspace2D_sptr
ProcessBackground::autoBackgroundSelection(Workspace2D_sptr bkgdWS) {
  // Get background type and create bakground function
  BackgroundFunction_sptr bkgdfunction = createBackgroundFunction(m_bkgdType);

  int bkgdorder = getProperty("BackgroundOrder");
  if (bkgdorder == 0)
    g_log.warning("(Input) background function order is 0.  It might not be "
                  "able to give a good estimation.");

  bkgdfunction->setAttributeValue("n", bkgdorder);
  bkgdfunction->initialize();

  g_log.information() << "Input background points has "
                      << bkgdWS->readX(0).size() << " data points for fit "
                      << bkgdorder << "-th order " << bkgdfunction->name()
                      << " (background) function" << bkgdfunction->asString()
                      << "\n";

  // Fit input (a few) background pionts to get initial guess
  API::IAlgorithm_sptr fit;
  try {
    fit = this->createChildAlgorithm("Fit", 0.0, 0.2, true);
  } catch (Exception::NotFoundError &) {
    g_log.error() << "Requires CurveFitting library." << std::endl;
    throw;
  }

  double startx = m_lowerBound;
  double endx = m_upperBound;
  fit->setProperty("Function",
                   boost::dynamic_pointer_cast<API::IFunction>(bkgdfunction));
  fit->setProperty("InputWorkspace", bkgdWS);
  fit->setProperty("WorkspaceIndex", 0);
  fit->setProperty("MaxIterations", 500);
  fit->setProperty("StartX", startx);
  fit->setProperty("EndX", endx);
  fit->setProperty("Minimizer", "Levenberg-Marquardt");
  fit->setProperty("CostFunction", "Least squares");

  fit->executeAsChildAlg();

  // Get fit result
  // a) Status
  std::string fitStatus = fit->getProperty("OutputStatus");
  bool allowedfailure = (fitStatus.find("cannot") < fitStatus.size()) &&
                        (fitStatus.find("tolerance") < fitStatus.size());
  if (fitStatus.compare("success") != 0 && !allowedfailure) {
    g_log.error() << "ProcessBackground: Fit Status = " << fitStatus
                  << ".  Not to update fit result" << std::endl;
    throw std::runtime_error("Bad Fit");
  }

  // b) check that chi2 got better
  const double chi2 = fit->getProperty("OutputChi2overDoF");
  g_log.information() << "Fit background: Fit Status = " << fitStatus
                      << ", chi2 = " << chi2 << "\n";

  // Filter and construct for the output workspace
  Workspace2D_sptr outws = filterForBackground(bkgdfunction);

  return outws;
} // END OF FUNCTION

/**
 * Creates and initialises a workspace with instrument definition and fills it
 * with data
 *
 * @param fileInfo information for the current file
 *
 * @param newFileNumber sequence number for the new file when added
 * into ws group
 *
 * @param buffer pre-allocated buffer to contain data values
 * @param imageY Object to set the Y data values in
 * @param imageE Object to set the E data values in
 *
 * @param parent A workspace which can be used to copy initialisation
 * information from (size/instrument def etc)
 *
 * @param loadAsRectImg if true, the new workspace will have one
 * spectrum per row and one bin per column, instead of the (default)
 * as many spectra as pixels.
 *
 * @param binSize size to rebin (1 == no re-bin == default)
 *
 * @param noiseThresh threshold for noise filtering
 *
 * @returns A newly created Workspace2D, as a shared pointer
 */
Workspace2D_sptr
LoadFITS::makeWorkspace(const FITSInfo &fileInfo, size_t &newFileNumber,
                        std::vector<char> &buffer, MantidImage &imageY,
                        MantidImage &imageE, const Workspace2D_sptr parent,
                        bool loadAsRectImg, int binSize, double noiseThresh) {
  // Create workspace (taking into account already here if rebinning is
  // going to happen)
  Workspace2D_sptr ws;
  if (!parent) {
    if (!loadAsRectImg) {
      size_t finalPixelCount = m_pixelCount / binSize * binSize;
      ws = boost::dynamic_pointer_cast<Workspace2D>(
          WorkspaceFactory::Instance().create("Workspace2D", finalPixelCount, 2,
                                              1));
    } else {
      ws = boost::dynamic_pointer_cast<Workspace2D>(
          WorkspaceFactory::Instance().create(
              "Workspace2D",
              fileInfo.axisPixelLengths[1] / binSize, // one bin per column
              fileInfo.axisPixelLengths[0] / binSize +
                  1, // one spectrum per row
              fileInfo.axisPixelLengths[0] / binSize));
    }
  } else {
    ws = boost::dynamic_pointer_cast<Workspace2D>(
        WorkspaceFactory::Instance().create(parent));
  }

  // this pixel scale property is used to set the workspace X values
  double cm_1 = getProperty("Scale");
  // amount of width units (e.g. cm) per pixel
  double cmpp = 1; // cm per pixel == bin width
  if (0.0 != cm_1)
    cmpp /= cm_1;
  cmpp *= static_cast<double>(binSize);

  if (loadAsRectImg && 1 == binSize) {
    // set data directly into workspace
    readDataToWorkspace(fileInfo, cmpp, ws, buffer);
  } else {
    readDataToImgs(fileInfo, imageY, imageE, buffer);
    doFilterNoise(noiseThresh, imageY, imageE);

    // Note this can change the sizes of the images and the number of pixels
    if (1 == binSize) {
      ws->setImageYAndE(imageY, imageE, 0, loadAsRectImg, cmpp,
                        false /* no parallel load */);

    } else {
      MantidImage rebinnedY(imageY.size() / binSize,
                            std::vector<double>(imageY[0].size() / binSize));
      MantidImage rebinnedE(imageE.size() / binSize,
                            std::vector<double>(imageE[0].size() / binSize));

      doRebin(binSize, imageY, imageE, rebinnedY, rebinnedE);
      ws->setImageYAndE(rebinnedY, rebinnedE, 0, loadAsRectImg, cmpp,
                        false /* no parallel load */);
    }
  }

  try {
    ws->setTitle(Poco::Path(fileInfo.filePath).getFileName());
  } catch (std::runtime_error &) {
    ws->setTitle(padZeros(newFileNumber, g_DIGIT_SIZE_APPEND));
  }
  ++newFileNumber;

  addAxesInfoAndLogs(ws, loadAsRectImg, fileInfo, binSize, cmpp);

  return ws;
}

  /** Construct output
   */
  Workspace2D_sptr RefinePowderInstrumentParameters2::genOutputWorkspace(FunctionDomain1DVector domain,
                                                               FunctionValues rawvalues)
  {
    // 1. Create and set up output workspace
    size_t lenx = m_dataWS->readX(m_wsIndex).size();
    size_t leny = m_dataWS->readY(m_wsIndex).size();

    Workspace2D_sptr outws = boost::dynamic_pointer_cast<Workspace2D>
        (WorkspaceFactory::Instance().create("Workspace2D", 6, lenx, leny));

    outws->getAxis(0)->setUnit("dSpacing");

    TextAxis* taxis = new TextAxis(outws->getNumberHistograms());
    taxis->setLabel(0, "Data");
    taxis->setLabel(1, "Model");
    taxis->setLabel(2, "DiffDM");
    taxis->setLabel(3, "Start");
    taxis->setLabel(4, "DiffDS");
    taxis->setLabel(5, "Zdiff");
    outws->replaceAxis(1, taxis);

    // 3. Re-calculate values
    FunctionValues funcvalues(domain);
    m_positionFunc->function(domain, funcvalues);

    // 4. Add values
    // a) X axis
    for (size_t iws = 0; iws < outws->getNumberHistograms(); ++iws)
    {
      MantidVec& vecX = outws->dataX(iws);
      for (size_t n = 0; n < lenx; ++n)
        vecX[n] = domain[n];
    }

    // b) Y axis
    const MantidVec& dataY = m_dataWS->readY(m_wsIndex);

    for (size_t i = 0; i < domain.size(); ++i)
    {
      outws->dataY(0)[i] = dataY[i];
      outws->dataY(1)[i] = funcvalues[i];
      outws->dataY(2)[i] = dataY[i] - funcvalues[i];
      outws->dataY(3)[i] = rawvalues[i];
      outws->dataY(4)[i] = dataY[i] - rawvalues[i];
    }

    // 5. Zscore
    vector<double> zscore = Kernel::getZscore(outws->readY(2));
    for (size_t i = 0; i < domain.size(); ++i)
      outws->dataY(5)[i] = zscore[i];

    return outws;
  }

  /** Executes the algorithm
   *
   *  @throw Exception::FileError If the grouping file cannot be opened or read successfully
   *  @throw runtime_error If unable to run one of the Child Algorithms successfully
   */
  void ReadGroupsFromFile::exec()
  {
    MatrixWorkspace_const_sptr ws = getProperty("InstrumentWorkspace");

    // Get the instrument.
    Instrument_const_sptr inst = ws->getInstrument();

    // Create a copy (without the data) of the workspace - it will contain the
    Workspace2D_sptr localWorkspace =
        boost::dynamic_pointer_cast<Workspace2D>(WorkspaceFactory::Instance().create(ws, ws->getNumberHistograms(), 2, 1));
    if (!localWorkspace)
      throw std::runtime_error("Failed when creating a Workspace2D from the input!");

    const std::string groupfile=getProperty("GroupingFilename");

    if ( ! groupfile.empty() )
    {
      std::string filename(groupfile);
      std::transform(filename.begin(), filename.end(), filename.begin(), tolower);
      if ( filename.find(".xml") != std::string::npos )
      {
        readXMLGroupingFile(groupfile);
      }
      else
      {
        readGroupingFile(groupfile);
      }
    }

    // Get the instrument.
    const int64_t nHist=localWorkspace->getNumberHistograms();

    // Determine whether the user wants to see unselected detectors or not
    const std::string su=getProperty("ShowUnselected");
    bool showunselected=(!su.compare("True"));
    bool success=false;

    for (int64_t i=0;i<nHist;i++)
    {
      ISpectrum * spec = localWorkspace->getSpectrum(i);
      const std::set<detid_t> & dets = spec->getDetectorIDs();
      if (dets.empty()) // Nothing
      {
        spec->dataY()[0]=0.0;
        continue;
      }
      // Find the first detector ID in the list
      calmap::const_iterator it=calibration.find(*dets.begin());
      if (it==calibration.end()) //Could not find the detector
      {
        spec->dataY()[0]=0.0;
        continue;
      }
      if (showunselected)
      {
        if (((*it).second).second==0)
          spec->dataY()[0]=0.0;
        else
          spec->dataY()[0]=static_cast<double>(((*it).second).first);
      }
      else
        spec->dataY()[0]=static_cast<double>(((*it).second).first);
      if (!success) success=true; //At least one detector is found in the cal file

    }
    progress(1);

    calibration.clear();
    if (!success) //Do some cleanup
    {
      localWorkspace.reset();
      throw std::runtime_error("Fail to found a detector in "+groupfile+" existing in instrument "+inst->getName());
    }
    setProperty("OutputWorkspace",localWorkspace);
    return;
  }

void TOFSANSResolution::exec()
{
  Workspace2D_sptr iqWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr reducedWS = getProperty("ReducedWorkspace");
  EventWorkspace_sptr reducedEventWS = boost::dynamic_pointer_cast<EventWorkspace>(reducedWS);
  const double min_wl = getProperty("MinWavelength");
  const double max_wl = getProperty("MaxWavelength");
  double pixel_size_x = getProperty("PixelSizeX");
  double pixel_size_y = getProperty("PixelSizeY");
  double R1 = getProperty("SourceApertureRadius");
  double R2 = getProperty("SampleApertureRadius");
  // Convert to meters
  pixel_size_x /= 1000.0;
  pixel_size_y /= 1000.0;
  R1 /= 1000.0;
  R2 /= 1000.0;
  wl_resolution = getProperty("DeltaT");

  // Although we want the 'ReducedWorkspace' to be an event workspace for this algorithm to do
  // anything, we don't want the algorithm to 'fail' if it isn't
  if (!reducedEventWS)
  {
    g_log.warning() << "An Event Workspace is needed to compute dQ. Calculation skipped." << std::endl;
    return;
  }

  // Calculate the output binning
  const std::vector<double> binParams = getProperty("OutputBinning");

  // Count histogram for normalization
  const int xLength = static_cast<int>(iqWS->readX(0).size());
  std::vector<double> XNorm(xLength-1, 0.0);

  // Create workspaces with each component of the resolution for debugging purposes
  MatrixWorkspace_sptr thetaWS = WorkspaceFactory::Instance().create(iqWS);
  declareProperty(new WorkspaceProperty<>("ThetaError","",Direction::Output));
  setPropertyValue("ThetaError","__"+iqWS->getName()+"_theta_error");
  setProperty("ThetaError",thetaWS);
  thetaWS->setX(0,iqWS->readX(0));
  MantidVec& ThetaY = thetaWS->dataY(0);

  MatrixWorkspace_sptr tofWS = WorkspaceFactory::Instance().create(iqWS);
  declareProperty(new WorkspaceProperty<>("TOFError","",Direction::Output));
  setPropertyValue("TOFError","__"+iqWS->getName()+"_tof_error");
  setProperty("TOFError",tofWS);
  tofWS->setX(0,iqWS->readX(0));
  MantidVec& TOFY = tofWS->dataY(0);

  // Initialize Dq
  MantidVec& DxOut = iqWS->dataDx(0);
  for ( int i = 0; i<xLength-1; i++ ) DxOut[i] = 0.0;

  const V3D samplePos = reducedWS->getInstrument()->getSample()->getPos();
  const V3D sourcePos = reducedWS->getInstrument()->getSource()->getPos();
  const V3D SSD = samplePos - sourcePos;
  const double L1 = SSD.norm();

  const int numberOfSpectra = static_cast<int>(reducedWS->getNumberHistograms());
  Progress progress(this,0.0,1.0,numberOfSpectra);

  PARALLEL_FOR2(reducedEventWS, iqWS)
  for (int i = 0; i < numberOfSpectra; i++)
  {
    PARALLEL_START_INTERUPT_REGION
    IDetector_const_sptr det;
    try {
      det = reducedEventWS->getDetector(i);
    } catch (Exception::NotFoundError&) {
      g_log.warning() << "Spectrum index " << i << " has no detector assigned to it - discarding" << std::endl;
      // 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 or if it's masked or a monitor, skip onto the next spectrum
    if ( !det || det->isMonitor() || det->isMasked() ) continue;

    // Get the flight path from the sample to the detector pixel
    const V3D scattered_flight_path = det->getPos() - samplePos;

    // Multiplicative factor to go from lambda to Q
    // Don't get fooled by the function name...
    const double theta = reducedEventWS->detectorTwoTheta(det);
    const double factor = 4.0 * M_PI * sin( theta/2.0 );

    EventList& el = reducedEventWS->getEventList(i);
    el.switchTo(WEIGHTED);

    std::vector<WeightedEvent>::iterator itev;
    std::vector<WeightedEvent>::iterator itev_end = el.getWeightedEvents().end();

    for (itev = el.getWeightedEvents().begin(); itev != itev_end; ++itev)
    {
      if ( itev->m_weight != itev->m_weight ) continue;
      if (std::abs(itev->m_weight) == std::numeric_limits<double>::infinity()) continue;
      if ( !isEmpty(min_wl) && itev->m_tof < min_wl ) continue;
      if ( !isEmpty(max_wl) && itev->m_tof > max_wl ) continue;

      const double q = factor/itev->m_tof;
      int iq = 0;

//.........这里部分代码省略.........

    /**
    * Executes the algorithm.
    * Saves the workspace specified by the user to the VTK XML format
    */
    void SaveVTK::exec()
    {
      std::string filename = getProperty("Filename");
      g_log.debug() << "Parameters: Filename='" << filename << "'" << std::endl;
      //add extension
      filename += ".vtu";

      MatrixWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
      if( !inputWorkspace )
      {
        g_log.error("Failed to retrieve inputWorkspace.");
        throw Exception::NullPointerException("SaveVTK::exec()", "inputWorkspace");
      }

      checkOptionalProperties();

      //Open file for writing
      std::ofstream outVTP(filename.c_str());
      if( !outVTP )
      {
        g_log.error("Failed to open file: " + filename);
        throw Exception::FileError("Failed to open file ", filename);
      }

      // First write document level XML header
      outVTP << "<?xml version=\"1.0\"?>\n"
        "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n"
        "<UnstructuredGrid>\n";

      const std::string workspaceID = inputWorkspace->id();    
      if( workspaceID.find("Workspace2D") != std::string::npos )
      {
        const Workspace2D_sptr localWorkspace = 
          boost::dynamic_pointer_cast<Workspace2D>(inputWorkspace);
        //      const size_t numberOfHist = localWorkspace->getNumberHistograms();

        //Write out whole range
        bool xMin(m_Xmin > 0.0), xMax(m_Xmax > 0.0);
        Progress prog(this,0.0,1.0,97);
        if( !xMin && !xMax )
        {
          for( int hNum = 2; hNum < 100; ++hNum )
          {
            writeVTKPiece(outVTP, localWorkspace->dataX(hNum), localWorkspace->dataY(hNum), 
              localWorkspace->dataE(hNum), hNum);
            prog.report();
          }
        }
        else
        {
          for( int hNum = 2; hNum < 100; ++hNum )
          {
            std::vector<double> xValue, yValue, errors;
            std::vector<double>::size_type nVals(localWorkspace->dataY(hNum).size());
            for( int i = 0; i < (int)nVals; ++i )
            {
              if( xMin && localWorkspace->dataX(hNum)[i] < m_Xmin ) continue;
              if( xMax && localWorkspace->dataX(hNum)[i+1] > m_Xmax)
              {
                xValue.push_back(localWorkspace->dataX(hNum)[i]);
                break;
              }
              xValue.push_back(localWorkspace->dataX(hNum)[i]);
              if( i == (int)nVals - 1 )
              {
                xValue.push_back(localWorkspace->dataX(hNum)[i+1]);
              } 	    
              yValue.push_back(localWorkspace->dataY(hNum)[i]);
              errors.push_back(localWorkspace->dataE(hNum)[i]);
            }
            //sanity check
            assert( (int)xValue.size() == (int)yValue.size() + 1 );

            writeVTKPiece(outVTP, xValue, yValue, errors, hNum);
            prog.report();
          }
        }
      }
      else
      {
        outVTP.close();
        Poco::File(filename).remove();
        throw Exception::NotImplementedError("SaveVTK only implemented for Workspace2D\n");
      }

      // Final XML end block tags
      outVTP << "</UnstructuredGrid>\n</VTKFile>\n";
      outVTP.close();
    }

void PoldiPeakSearch::exec() {
  g_log.information() << "PoldiPeakSearch:" << std::endl;

  Workspace2D_sptr correlationWorkspace = getProperty("InputWorkspace");
  MantidVec correlationQValues = correlationWorkspace->readX(0);
  MantidVec correlatedCounts = correlationWorkspace->readY(0);
  g_log.information() << "   Auto-correlation data read." << std::endl;

  Unit_sptr xUnit = correlationWorkspace->getAxis(0)->unit();

  if (xUnit->caption() == "") {
    g_log.information()
        << "   Workspace does not have unit, defaulting to MomentumTransfer."
        << std::endl;

    xUnit = UnitFactory::Instance().create("MomentumTransfer");
  } else {
    g_log.information() << "   Unit of workspace is " << xUnit->caption() << "."
                        << std::endl;
  }

  setMinimumDistance(getProperty("MinimumPeakSeparation"));
  setMinimumPeakHeight(getProperty("MinimumPeakHeight"));
  setMaximumPeakNumber(getProperty("MaximumPeakNumber"));

  if (m_doubleMinimumDistance > static_cast<int>(correlatedCounts.size())) {
    throw(std::runtime_error("MinimumPeakSeparation is smaller than number of "
                             "spectrum points - no peaks possible."));
  }

  g_log.information() << "   Parameters set." << std::endl;

  MantidVec summedNeighborCounts = getNeighborSums(correlatedCounts);
  g_log.information() << "   Neighboring counts summed, contains "
                      << summedNeighborCounts.size() << " data points."
                      << std::endl;

  std::list<MantidVec::const_iterator> peakPositionsSummed =
      findPeaks(summedNeighborCounts.begin(), summedNeighborCounts.end());
  g_log.information() << "   Peaks detected in summed spectrum: "
                      << peakPositionsSummed.size() << std::endl;

  /* This step is required because peaks are actually searched in the
   * "sum-of-neighbors"-spectrum.
   * The mapping removes the offset from the peak position which results from
   * different beginning
   * of this vector compared to the original correlation counts.
   */
  std::list<MantidVec::const_iterator> peakPositionsCorrelation =
      mapPeakPositionsToCorrelationData(peakPositionsSummed,
                                        summedNeighborCounts.begin(),
                                        correlatedCounts.begin());
  g_log.information() << "   Peak positions transformed to original spectrum."
                      << std::endl;

  /* Since intensities are required for filtering, they are extracted from the
   * original count data,
   * along with the Q-values.
   */
  std::vector<PoldiPeak_sptr> peakCoordinates =
      getPeaks(correlatedCounts.begin(), correlatedCounts.end(),
               peakPositionsCorrelation, correlationQValues, xUnit);
  g_log.information()
      << "   Extracted peak positions in Q and intensity guesses." << std::endl;

  UncertainValue backgroundWithSigma =
      getBackgroundWithSigma(peakPositionsCorrelation, correlatedCounts);
  g_log.information() << "   Calculated average background and deviation: "
                      << UncertainValueIO::toString(backgroundWithSigma)
                      << std::endl;

  if ((*getProperty("MinimumPeakHeight")).isDefault()) {
    setMinimumPeakHeight(minimumPeakHeightFromBackground(backgroundWithSigma));
  }

  std::vector<PoldiPeak_sptr> intensityFilteredPeaks(peakCoordinates.size());
  auto newEnd = std::remove_copy_if(
      peakCoordinates.begin(), peakCoordinates.end(),
      intensityFilteredPeaks.begin(),
      boost::bind(&PoldiPeakSearch::isLessThanMinimum, this, _1));
  intensityFilteredPeaks.resize(
      std::distance(intensityFilteredPeaks.begin(), newEnd));

  g_log.information() << "   Peaks above minimum intensity ("
                      << m_minimumPeakHeight
                      << "): " << intensityFilteredPeaks.size() << std::endl;

  std::sort(intensityFilteredPeaks.begin(), intensityFilteredPeaks.end(),
            boost::bind<bool>(&PoldiPeak::greaterThan, _1, _2,
                              &PoldiPeak::intensity));

  for (std::vector<PoldiPeak_sptr>::const_iterator peak =
           intensityFilteredPeaks.begin();
       peak != intensityFilteredPeaks.end(); ++peak) {
    m_peaks->addPeak(*peak);
  }

  /* The derived background error is set as error in the workspace containing
   * correlation data, so it may be used as weights for peak fitting later on.
   */
//.........这里部分代码省略.........