IntegratePeaksMD2.cpp

#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidMDAlgorithms/GSLFunctions.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidDataObjects/PeakShapeSpherical.h"
#include "MantidKernel/System.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidMDAlgorithms/IntegratePeaksMD2.h"
#include "MantidDataObjects/CoordTransformDistance.h"
#include "MantidKernel/ListValidator.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/TextAxis.h"
#include "MantidKernel/Utils.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/Column.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/Progress.h"
#include <boost/math/special_functions/fpclassify.hpp>
#include <gsl/gsl_integration.h>
#include <fstream>

namespace Mantid {
namespace MDAlgorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IntegratePeaksMD2)

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;

//----------------------------------------------------------------------------------------------
/** Constructor
 */
IntegratePeaksMD2::IntegratePeaksMD2() {}

//----------------------------------------------------------------------------------------------
/** Destructor
 */
IntegratePeaksMD2::~IntegratePeaksMD2() {}

//----------------------------------------------------------------------------------------------

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void IntegratePeaksMD2::init() {
  declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace", "",
                                                           Direction::Input),
                  "An input MDEventWorkspace.");

  std::vector<std::string> propOptions;
  propOptions.push_back("Q (lab frame)");
  propOptions.push_back("Q (sample frame)");
  propOptions.push_back("HKL");

  declareProperty(
      new PropertyWithValue<double>("PeakRadius", 1.0, Direction::Input),
      "Fixed radius around each peak position in which to integrate (in the "
      "same units as the workspace).");

  declareProperty(
      new PropertyWithValue<double>("BackgroundInnerRadius", 0.0,
                                    Direction::Input),
      "Inner radius to use to evaluate the background of the peak.\n"
      "If smaller than PeakRadius, then we assume BackgroundInnerRadius = "
      "PeakRadius.");

  declareProperty(
      new PropertyWithValue<double>("BackgroundOuterRadius", 0.0,
                                    Direction::Input),
      "Outer radius to use to evaluate the background of the peak.\n"
      "The signal density around the peak (BackgroundInnerRadius < r < "
      "BackgroundOuterRadius) is used to estimate the background under the "
      "peak.\n"
      "If smaller than PeakRadius, no background measurement is done.");

  declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace", "",
                                                        Direction::Input),
                  "A PeaksWorkspace containing the peaks to integrate.");

  declareProperty(
      new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace", "",
                                            Direction::Output),
      "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
      "with the peaks' integrated intensities.");

  declareProperty("ReplaceIntensity", true,
                  "Always replace intensity in PeaksWorkspacem (default).\n"
                  "If false, then do not replace intensity if calculated value "
                  "is 0 (used for SNSSingleCrystalReduction)");

  declareProperty(
      "IntegrateIfOnEdge", true,
      "Only warning if all of peak outer radius is not on detector (default).\n"
      "If false, do not integrate if the outer radius is not on a detector.");

  declareProperty("AdaptiveQBackground", false,
                  "Default is false.   If true, all background values"
                  "vary on a line so that they are"
                  "background plus AdaptiveQMultiplier multiplied"
                  "by the magnitude of Q at the peak center so each peak has a "
                  "different integration radius.  Q includes the 2*pi factor.");

  declareProperty("Cylinder", false,
                  "Default is sphere.  Use next five parameters for cylinder.");

  declareProperty(
      new PropertyWithValue<double>("CylinderLength", 0.0, Direction::Input),
      "Length of cylinder in which to integrate (in the same units as the "
      "workspace).");

  declareProperty(
      new PropertyWithValue<double>("PercentBackground", 0.0, Direction::Input),
      "Percent of CylinderLength that is background (20 is 20%)");

  std::vector<std::string> peakNames =
      FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
  peakNames.push_back("NoFit");
  declareProperty("ProfileFunction", "Gaussian",
                  boost::make_shared<StringListValidator>(peakNames),
                  "Fitting function for profile that is used only with "
                  "Cylinder integration.");

  std::vector<std::string> integrationOptions(2);
  integrationOptions[0] = "Sum";
  integrationOptions[1] = "GaussianQuadrature";
  auto integrationvalidator =
      boost::make_shared<StringListValidator>(integrationOptions);
  declareProperty("IntegrationOption", "GaussianQuadrature",
                  integrationvalidator,
                  "Integration method for calculating intensity "
                  "used only with Cylinder integration.");

  declareProperty(
      new FileProperty("ProfilesFile", "", FileProperty::OptionalSave,
                       std::vector<std::string>(1, "profiles")),
      "Save (Optionally) as Isaw peaks file with profiles included");

  declareProperty("AdaptiveQMultiplier", 0.0,
                  "Peak integration radius varies on a line so that it is"
                  "PeakRadius plus this value multiplied"
                  "by the magnitude of Q at the peak center so each peak has a "
                  "different integration radius.  Q includes the 2*pi factor.");
}

//----------------------------------------------------------------------------------------------
/** Integrate the peaks of the workspace using parameters saved in the algorithm
 * class
 * @param ws ::  MDEventWorkspace to integrate
 */
template <typename MDE, size_t nd>
void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
  if (nd != 3)
    throw std::invalid_argument("For now, we expect the input MDEventWorkspace "
                                "to have 3 dimensions only.");

  /// Peak workspace to integrate
  Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS =
      getProperty("PeaksWorkspace");

  /// Output peaks workspace, create if needed
  Mantid::DataObjects::PeaksWorkspace_sptr peakWS =
      getProperty("OutputWorkspace");
  if (peakWS != inPeakWS)
    peakWS.reset(inPeakWS->clone().release());
  // This only fails in the unit tests which say that MaskBTP is not registered
  try {
    runMaskDetectors(inPeakWS, "Tube", "edges");
    runMaskDetectors(inPeakWS, "Pixel", "edges");
  } catch (...) {
    g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
                "edge for IntegrateIfOnEdge option");
  }

  // Get the instrument and its detectors
  Geometry::Instrument_const_sptr inst = inPeakWS->getInstrument();
  calculateE1(inst); // fill E1Vec for use in detectorQ
  Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse =
      ws->getSpecialCoordinateSystem();

  /// Radius to use around peaks
  double PeakRadius = getProperty("PeakRadius");
  /// Background (end) radius
  double BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
  /// Start radius of the background
  double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
  if (BackgroundInnerRadius < PeakRadius)
    BackgroundInnerRadius = PeakRadius;
  /// Cylinder Length to use around peaks for cylinder
  double cylinderLength = getProperty("CylinderLength");
  Workspace2D_sptr wsProfile2D, wsFit2D, wsDiff2D;
  size_t numSteps = 0;
  bool cylinderBool = getProperty("Cylinder");
  bool adaptiveQBackground = getProperty("AdaptiveQBackground");
  double adaptiveQMultiplier = getProperty("AdaptiveQMultiplier");
  double adaptiveQBackgroundMultiplier = 0.0;
  if (adaptiveQBackground)
    adaptiveQBackgroundMultiplier = adaptiveQMultiplier;
  std::vector<double> PeakRadiusVector(peakWS->getNumberPeaks(), PeakRadius);
  std::vector<double> BackgroundInnerRadiusVector(peakWS->getNumberPeaks(),
                                                  BackgroundInnerRadius);
  std::vector<double> BackgroundOuterRadiusVector(peakWS->getNumberPeaks(),
                                                  BackgroundOuterRadius);
  if (cylinderBool) {
    numSteps = 100;
    size_t histogramNumber = peakWS->getNumberPeaks();
    Workspace_sptr wsProfile = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
    wsProfile2D = boost::dynamic_pointer_cast<Workspace2D>(wsProfile);
    AnalysisDataService::Instance().addOrReplace("ProfilesData", wsProfile2D);
    Workspace_sptr wsFit = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
    wsFit2D = boost::dynamic_pointer_cast<Workspace2D>(wsFit);
    AnalysisDataService::Instance().addOrReplace("ProfilesFit", wsFit2D);
    Workspace_sptr wsDiff = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
    wsDiff2D = boost::dynamic_pointer_cast<Workspace2D>(wsDiff);
    AnalysisDataService::Instance().addOrReplace("ProfilesFitDiff", wsDiff2D);
    TextAxis *const newAxis1 = new TextAxis(peakWS->getNumberPeaks());
    TextAxis *const newAxis2 = new TextAxis(peakWS->getNumberPeaks());
    TextAxis *const newAxis3 = new TextAxis(peakWS->getNumberPeaks());
    wsProfile2D->replaceAxis(1, newAxis1);
    wsFit2D->replaceAxis(1, newAxis2);
    wsDiff2D->replaceAxis(1, newAxis3);
    for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
      // Get a direct ref to that peak.
      IPeak &p = peakWS->getPeak(i);
      std::ostringstream label;
      label << Utils::round(p.getH()) << "_" << Utils::round(p.getK()) << "_"
            << Utils::round(p.getL()) << "_" << p.getRunNumber();
      newAxis1->setLabel(i, label.str());
      newAxis2->setLabel(i, label.str());
      newAxis3->setLabel(i, label.str());
    }
  }
  double percentBackground = getProperty("PercentBackground");
  size_t peakMin = 0;
  size_t peakMax = numSteps;
  double ratio = 0.0;
  if (cylinderBool) {
    peakMin = static_cast<size_t>(static_cast<double>(numSteps) *
                                  percentBackground / 100.);
    peakMax = numSteps - peakMin - 1;
    size_t numPeakCh = peakMax - peakMin + 1; // number of peak channels
    size_t numBkgCh = numSteps - numPeakCh;   // number of background channels
    ratio = static_cast<double>(numPeakCh) / static_cast<double>(numBkgCh);
  }
  /// Replace intensity with 0
  bool replaceIntensity = getProperty("ReplaceIntensity");
  bool integrateEdge = getProperty("IntegrateIfOnEdge");

  std::string profileFunction = getProperty("ProfileFunction");
  std::string integrationOption = getProperty("IntegrationOption");
  std::ofstream out;
  if (cylinderBool && profileFunction.compare("NoFit") != 0) {
    std::string outFile = getProperty("InputWorkspace");
    outFile.append(profileFunction);
    outFile.append(".dat");
    std::string save_path =
        ConfigService::Instance().getString("defaultsave.directory");
    outFile = save_path + outFile;
    out.open(outFile.c_str(), std::ofstream::out);
  }
  //
  // If the following OMP pragma is included, this algorithm seg faults
  // sporadically when processing multiple TOPAZ runs in a script, on
  // Scientific Linux 6.2.  Typically, it seg faults after 2 to 6 runs are
  // processed, though occasionally it will process all 8 requested in the
  // script without crashing.  Since the lower level codes already use OpenMP,
  // parallelizing at this level is only marginally useful, giving about a
  // 5-10% speedup.  Perhaps is should just be removed permanantly, but for
  // now it is commented out to avoid the seg faults.  Refs #5533
  // PRAGMA_OMP(parallel for schedule(dynamic, 10) )
  // Initialize progress reporting
  int nPeaks = peakWS->getNumberPeaks();
  Progress progress(this, 0., 1., nPeaks);
  for (int i = 0; i < nPeaks; ++i) {
    if (this->getCancel())
      break; // User cancellation
    progress.report();

    // Get a direct ref to that peak.
    IPeak &p = peakWS->getPeak(i);

    // Get the peak center as a position in the dimensions of the workspace
    V3D pos;
    if (CoordinatesToUse == Mantid::Kernel::QLab) //"Q (lab frame)"
      pos = p.getQLabFrame();
    else if (CoordinatesToUse == Mantid::Kernel::QSample) //"Q (sample frame)"
      pos = p.getQSampleFrame();
    else if (CoordinatesToUse == Mantid::Kernel::HKL) //"HKL"
      pos = p.getHKL();

    // Do not integrate if sphere is off edge of detector

    if (!detectorQ(p.getQLabFrame(),
                   std::max(BackgroundOuterRadius, PeakRadius))) {
      g_log.warning() << "Warning: sphere/cylinder for integration is off edge "
                         "of detector for peak " << i << std::endl;
      if (!integrateEdge) {
        if (replaceIntensity) {
          p.setIntensity(0.0);
          p.setSigmaIntensity(0.0);
        }
        continue;
      }
    }

    // Build the sphere transformation
    bool dimensionsUsed[nd];
    coord_t center[nd];
    for (size_t d = 0; d < nd; ++d) {
      dimensionsUsed[d] = true; // Use all dimensions
      center[d] = static_cast<coord_t>(pos[d]);
    }
    signal_t signal = 0;
    signal_t errorSquared = 0;
    signal_t bgSignal = 0;
    signal_t bgErrorSquared = 0;
    double background_total = 0.0;
    if (!cylinderBool) {
      // modulus of Q
      coord_t lenQpeak = 0.0;
      if (adaptiveQMultiplier > 0.0) {
        lenQpeak = 0.0;
        for (size_t d = 0; d < nd; d++) {
          lenQpeak += center[d] * center[d];
        }
        lenQpeak = std::sqrt(lenQpeak);
      }

      PeakRadiusVector[i] = adaptiveQMultiplier * lenQpeak + PeakRadius;
      BackgroundInnerRadiusVector[i] =
          adaptiveQBackgroundMultiplier * lenQpeak + BackgroundInnerRadius;
      BackgroundOuterRadiusVector[i] =
          adaptiveQBackgroundMultiplier * lenQpeak + BackgroundOuterRadius;
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      if (Peak *shapeablePeak = dynamic_cast<Peak *>(&p)) {

        PeakShape *sphere = new PeakShapeSpherical(
            PeakRadiusVector[i], BackgroundInnerRadiusVector[i],
            BackgroundOuterRadiusVector[i], CoordinatesToUse, this->name(),
            this->version());
        shapeablePeak->setPeakShape(sphere);
      }

      // Perform the integration into whatever box is contained within.
      ws->getBox()->integrateSphere(
          sphere,
          static_cast<coord_t>((adaptiveQMultiplier * lenQpeak + PeakRadius) *
                               (adaptiveQMultiplier * lenQpeak + PeakRadius)),
          signal, errorSquared);

      // Integrate around the background radius

      if (BackgroundOuterRadius > PeakRadius) {
        // Get the total signal inside "BackgroundOuterRadius"
        ws->getBox()->integrateSphere(
            sphere,
            static_cast<coord_t>((adaptiveQBackgroundMultiplier * lenQpeak +
                                  BackgroundOuterRadius) *
                                 (adaptiveQBackgroundMultiplier * lenQpeak +
                                  BackgroundOuterRadius)),
            bgSignal, bgErrorSquared);

        // Evaluate the signal inside "BackgroundInnerRadius"
        signal_t interiorSignal = 0;
        signal_t interiorErrorSquared = 0;

        // Integrate this 3rd radius, if needed
        if (BackgroundInnerRadius != PeakRadius) {
          ws->getBox()->integrateSphere(
              sphere,
              static_cast<coord_t>((adaptiveQBackgroundMultiplier * lenQpeak +
                                    BackgroundInnerRadius) *
                                   (adaptiveQBackgroundMultiplier * lenQpeak +
                                    BackgroundInnerRadius)),
              interiorSignal, interiorErrorSquared);
        } else {
          // PeakRadius == BackgroundInnerRadius, so use the previous value
          interiorSignal = signal;
          interiorErrorSquared = errorSquared;
        }
        // Subtract the peak part to get the intensity in the shell
        // (BackgroundInnerRadius < r < BackgroundOuterRadius)
        bgSignal -= interiorSignal;
        // We can subtract the error (instead of adding) because the two values
        // are 100% dependent; this is the same as integrating a shell.
        bgErrorSquared -= interiorErrorSquared;

        // Relative volume of peak vs the BackgroundOuterRadius sphere
        double ratio = (PeakRadius / BackgroundOuterRadius);
        double peakVolume = ratio * ratio * ratio;

        // Relative volume of the interior of the shell vs overall background
        double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
        // Volume of the bg shell, relative to the volume of the
        // BackgroundOuterRadius sphere
        double bgVolume = 1.0 - interiorRatio * interiorRatio * interiorRatio;

        // Finally, you will multiply the bg intensity by this to get the
        // estimated background under the peak volume
        double scaleFactor = peakVolume / bgVolume;
        bgSignal *= scaleFactor;
        bgErrorSquared *= scaleFactor * scaleFactor;
      }
    } else {
      CoordTransformDistance cylinder(nd, center, dimensionsUsed, 2);

      // Perform the integration into whatever box is contained within.
      std::vector<signal_t> signal_fit;

      signal_fit.clear();
      for (size_t j = 0; j < numSteps; j++)
        signal_fit.push_back(0.0);
      ws->getBox()->integrateCylinder(cylinder,
                                      static_cast<coord_t>(PeakRadius),
                                      static_cast<coord_t>(cylinderLength),
                                      signal, errorSquared, signal_fit);
      for (size_t j = 0; j < numSteps; j++) {
        wsProfile2D->dataX(i)[j] = static_cast<double>(j);
        wsProfile2D->dataY(i)[j] = signal_fit[j];
        wsProfile2D->dataE(i)[j] = std::sqrt(signal_fit[j]);
      }

      // Integrate around the background radius
      if (BackgroundOuterRadius > PeakRadius) {
        // Get the total signal inside "BackgroundOuterRadius"

        signal_fit.clear();
        for (size_t j = 0; j < numSteps; j++)
          signal_fit.push_back(0.0);
        ws->getBox()->integrateCylinder(
            cylinder, static_cast<coord_t>(BackgroundOuterRadius),
            static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
            signal_fit);
        for (size_t j = 0; j < numSteps; j++) {
          wsProfile2D->dataX(i)[j] = static_cast<double>(j);
          wsProfile2D->dataY(i)[j] = signal_fit[j];
          wsProfile2D->dataE(i)[j] = std::sqrt(signal_fit[j]);
        }

        // Evaluate the signal inside "BackgroundInnerRadius"
        signal_t interiorSignal = 0;
        signal_t interiorErrorSquared = 0;

        // Integrate this 3rd radius, if needed
        if (BackgroundInnerRadius != PeakRadius) {
          ws->getBox()->integrateCylinder(
              cylinder, static_cast<coord_t>(BackgroundInnerRadius),
              static_cast<coord_t>(cylinderLength), interiorSignal,
              interiorErrorSquared, signal_fit);
        } else {
          // PeakRadius == BackgroundInnerRadius, so use the previous value
          interiorSignal = signal;
          interiorErrorSquared = errorSquared;
        }
        // Subtract the peak part to get the intensity in the shell
        // (BackgroundInnerRadius < r < BackgroundOuterRadius)
        bgSignal -= interiorSignal;
        // We can subtract the error (instead of adding) because the two values
        // are 100% dependent; this is the same as integrating a shell.
        bgErrorSquared -= interiorErrorSquared;
        // Relative volume of peak vs the BackgroundOuterRadius cylinder
        double ratio = (PeakRadius / BackgroundOuterRadius);
        double peakVolume = ratio * ratio * cylinderLength;

        // Relative volume of the interior of the shell vs overall background
        double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
        // Volume of the bg shell, relative to the volume of the
        // BackgroundOuterRadius cylinder
        double bgVolume = 1.0 - interiorRatio * interiorRatio * cylinderLength;

        // Finally, you will multiply the bg intensity by this to get the
        // estimated background under the peak volume
        double scaleFactor = peakVolume / bgVolume;
        bgSignal *= scaleFactor;
        bgErrorSquared *= scaleFactor * scaleFactor;
      } else {
        for (size_t j = 0; j < numSteps; j++) {
          wsProfile2D->dataX(i)[j] = static_cast<double>(j);
          wsProfile2D->dataY(i)[j] = signal_fit[j];
          wsProfile2D->dataE(i)[j] = std::sqrt(signal_fit[j]);
        }
      }

      if (profileFunction.compare("NoFit") == 0) {
        signal = 0.;
        for (size_t j = 0; j < numSteps; j++) {
          if (j < peakMin || j > peakMax)
            background_total = background_total + wsProfile2D->dataY(i)[j];
          else
            signal = signal + wsProfile2D->dataY(i)[j];
        }
        errorSquared = std::fabs(signal);
      } else {
        API::IAlgorithm_sptr findpeaks =
            createChildAlgorithm("FindPeaks", -1, -1, false);
        findpeaks->setProperty("InputWorkspace", wsProfile2D);
        findpeaks->setProperty<int>("FWHM", 7);
        findpeaks->setProperty<int>("Tolerance", 4);
        // FindPeaks will do the checking on the validity of WorkspaceIndex
        findpeaks->setProperty("WorkspaceIndex", static_cast<int>(i));

        // Get the specified peak positions, which is optional
        findpeaks->setProperty<std::string>("PeakFunction", profileFunction);
        // FindPeaks will use linear or flat if they are better
        findpeaks->setProperty<std::string>("BackgroundType", "Quadratic");
        findpeaks->setProperty<bool>("HighBackground", true);
        findpeaks->setProperty<bool>("RawPeakParameters", true);
        std::vector<double> peakPosToFit;
        peakPosToFit.push_back(static_cast<double>(numSteps) / 2.0);
        findpeaks->setProperty("PeakPositions", peakPosToFit);
        findpeaks->setProperty<int>("MinGuessedPeakWidth", 4);
        findpeaks->setProperty<int>("MaxGuessedPeakWidth", 4);
        try {
          findpeaks->executeAsChildAlg();
        } catch (...) {
          g_log.error("Can't execute FindPeaks algorithm");
          continue;
        }

        API::ITableWorkspace_sptr paramws = findpeaks->getProperty("PeaksList");
        if (paramws->rowCount() < 1)
          continue;
        std::ostringstream fun_str;
        fun_str << "name=" << profileFunction;

        size_t numcols = paramws->columnCount();
        std::vector<std::string> paramsName = paramws->getColumnNames();
        std::vector<double> paramsValue;
        API::TableRow row = paramws->getRow(0);
        int spectrum;
        row >> spectrum;
        for (size_t j = 1; j < numcols; ++j) {
          double parvalue;
          row >> parvalue;
          if (j == numcols - 4)
            fun_str << ";name=Quadratic";
          // erase f0. or f1.
          // if (j > 0 && j < numcols-1) fun_str << "," <<
          // paramsName[j].erase(0,3) <<"="<<parvalue;
          if (j > 0 && j < numcols - 1)
            fun_str << "," << paramsName[j] << "=" << parvalue;
          paramsValue.push_back(parvalue);
        }
        if (i == 0) {
          for (size_t j = 0; j < numcols; ++j)
            out << std::setw(20) << paramsName[j] << " ";
          out << "\n";
        }
        out << std::setw(20) << i;
        for (size_t j = 0; j < numcols - 1; ++j)
          out << std::setw(20) << std::fixed << std::setprecision(10)
              << paramsValue[j] << " ";
        out << "\n";

        // Evaluate fit at points

        IFunction_sptr ifun =
            FunctionFactory::Instance().createInitialized(fun_str.str());
        boost::shared_ptr<const CompositeFunction> fun =
            boost::dynamic_pointer_cast<const CompositeFunction>(ifun);
        const Mantid::MantidVec &x = wsProfile2D->readX(i);
        wsFit2D->dataX(i) = x;
        wsDiff2D->dataX(i) = x;
        FunctionDomain1DVector domain(x);
        FunctionValues yy(domain);
        fun->function(domain, yy);
        const Mantid::MantidVec &yValues = wsProfile2D->readY(i);
        for (size_t j = 0; j < numSteps; j++) {
          wsFit2D->dataY(i)[j] = yy[j];
          wsDiff2D->dataY(i)[j] = yValues[j] - yy[j];
        }

        // Calculate intensity
        signal = 0.0;
        if (integrationOption.compare("Sum") == 0) {
          for (size_t j = peakMin; j <= peakMax; j++)
            if (!boost::math::isnan(yy[j]) && !boost::math::isinf(yy[j]))
              signal += yy[j];
        } else {
          gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);

          double error;

          gsl_function F;
          F.function = &Mantid::MDAlgorithms::f_eval2;
          F.params = &fun;

          gsl_integration_qags(&F, x[peakMin], x[peakMax], 0, 1e-7, 1000, w,
                               &signal, &error);

          gsl_integration_workspace_free(w);
        }
        errorSquared = std::fabs(signal);
        // Get background counts
        for (size_t j = 0; j < numSteps; j++) {
          // paramsValue[numcols-2] is chisq
          double background = paramsValue[numcols - 3] * x[j] * x[j] +
                              paramsValue[numcols - 4] * x[j] +
                              paramsValue[numcols - 5];
          if (j < peakMin || j > peakMax)
            background_total = background_total + background;
        }
      }
    }
    checkOverlap(
        i, peakWS, CoordinatesToUse,
        2.0 * std::max(PeakRadiusVector[i], BackgroundOuterRadiusVector[i]));
    // Save it back in the peak object.
    if (signal != 0. || replaceIntensity) {
      p.setIntensity(signal - ratio * background_total - bgSignal);
      p.setSigmaIntensity(sqrt(errorSquared +
                               ratio * ratio * std::fabs(background_total) +
                               bgErrorSquared));
    }

    g_log.information() << "Peak " << i << " at " << pos << ": signal "
                        << signal << " (sig^2 " << errorSquared
                        << "), with background "
                        << bgSignal + ratio * background_total << " (sig^2 "
                        << bgErrorSquared +
                               ratio * ratio * std::fabs(background_total)
                        << ") subtracted." << std::endl;
  }
  // This flag is used by the PeaksWorkspace to evaluate whether it has been
  // integrated.
  peakWS->mutableRun().addProperty("PeaksIntegrated", 1, true);
  // These flags are specific to the algorithm.
  peakWS->mutableRun().addProperty("PeakRadius", PeakRadiusVector, true);
  peakWS->mutableRun().addProperty("BackgroundInnerRadius",
                                   BackgroundInnerRadiusVector, true);
  peakWS->mutableRun().addProperty("BackgroundOuterRadius",
                                   BackgroundOuterRadiusVector, true);

  // save profiles in peaks file
  const std::string outfile = getProperty("ProfilesFile");
  if (outfile.length() > 0) {
    IAlgorithm_sptr alg;
    try {
      alg = createChildAlgorithm("SaveIsawPeaks", -1, -1, false);
    } catch (Exception::NotFoundError &) {
      g_log.error("Can't locate SaveIsawPeaks algorithm");
      throw;
    }
    alg->setProperty("InputWorkspace", peakWS);
    alg->setProperty("ProfileWorkspace", wsProfile2D);
    alg->setPropertyValue("Filename", outfile);
    alg->execute();
  }
  // Save the output
  setProperty("OutputWorkspace", peakWS);
}

/*
 * Define edges for each instrument by masking. For CORELLI, tubes 1 and 16, and
 *pixels 0 and 255.
 * Get Q in the lab frame for every peak, call it C
 * For every point on the edge, the trajectory in reciprocal space is a straight
 *line, going through O=V3D(0,0,0).
 * Calculate a point at a fixed momentum, say k=1. Q in the lab frame
 *E=V3D(-k*sin(tt)*cos(ph),-k*sin(tt)*sin(ph),k-k*cos(ph)).
 * Normalize E to 1: E=E*(1./E.norm())
 *
 * @param inst: instrument
 */
void IntegratePeaksMD2::calculateE1(Geometry::Instrument_const_sptr inst) {
  std::vector<detid_t> detectorIDs = inst->getDetectorIDs();

  for (auto detID = detectorIDs.begin(); detID != detectorIDs.end(); ++detID) {
    Mantid::Geometry::IDetector_const_sptr det = inst->getDetector(*detID);
    if (det->isMonitor())
      continue; // skip monitor
    if (!det->isMasked())
      continue; // edge is masked so don't check if not masked
    double tt1 = det->getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)); // two theta
    double ph1 = det->getPhi();                                // phi
    V3D E1 = V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
                 1. - std::cos(tt1)); // end of trajectory
    E1 = E1 * (1. / E1.norm());       // normalize
    E1Vec.push_back(E1);
  }
}

/** Calculate if this Q is on a detector
 * The distance from C to OE is given by dv=C-E*(C.scalar_prod(E))
 * If dv.norm<integration_radius, one of the detector trajectories on the edge
 *is too close to the peak
 * This method is applied to all masked pixels. If there are masked pixels
 *trajectories inside an integration volume, the peak must be rejected.
 *
 * @param QLabFrame: The Peak center.
 * @param r: Peak radius.
 */
bool IntegratePeaksMD2::detectorQ(Mantid::Kernel::V3D QLabFrame, double r) {

  for (auto E1 = E1Vec.begin(); E1 != E1Vec.end(); ++E1) {
    V3D distv = QLabFrame -
                *E1 * (QLabFrame.scalar_prod(
                          *E1)); // distance to the trajectory as a vector
    if (distv.norm() < r) {
      return false;
    }
  }
  return true;
}
void IntegratePeaksMD2::runMaskDetectors(
    Mantid::DataObjects::PeaksWorkspace_sptr peakWS, std::string property,
    std::string values) {
  IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
  alg->setProperty<Workspace_sptr>("Workspace", peakWS);
  alg->setProperty(property, values);
  if (!alg->execute())
    throw std::runtime_error(
        "MaskDetectors Child Algorithm has not executed successfully");
}

void IntegratePeaksMD2::checkOverlap(
    int i, Mantid::DataObjects::PeaksWorkspace_sptr peakWS,
    Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse, double radius) {
  // Get a direct ref to that peak.
  IPeak &p1 = peakWS->getPeak(i);
  V3D pos1;
  if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
    pos1 = p1.getQLabFrame();
  else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
    pos1 = p1.getQSampleFrame();
  else if (CoordinatesToUse == Kernel::HKL) //"HKL"
    pos1 = p1.getHKL();
  for (int j = i + 1; j < peakWS->getNumberPeaks(); ++j) {
    // Get a direct ref to rest of peaks peak.
    IPeak &p2 = peakWS->getPeak(j);
    V3D pos2;
    if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
      pos2 = p2.getQLabFrame();
    else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
      pos2 = p2.getQSampleFrame();
    else if (CoordinatesToUse == Kernel::HKL) //"HKL"
      pos2 = p2.getHKL();
    if (pos1.distance(pos2) < radius) {
      g_log.warning() << " Warning:  Peak integration spheres for peaks " << i
                      << " and " << j << " overlap.  Distance between peaks is "
                      << pos1.distance(pos2) << std::endl;
    }
  }
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void IntegratePeaksMD2::exec() {
  inWS = getProperty("InputWorkspace");

  CALL_MDEVENT_FUNCTION(this->integrate, inWS);
}

double f_eval2(double x, void *params) {
  boost::shared_ptr<const API::CompositeFunction> fun =
      *(boost::shared_ptr<const API::CompositeFunction> *)params;
  FunctionDomain1DVector domain(x);
  FunctionValues yval(domain);
  fun->function(domain, yval);
  return yval[0];
}

} // namespace Mantid
} // namespace DataObjects