//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/CalculateEfficiency.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/EventList.h"
#include <vector>
#include "MantidKernel/BoundedValidator.h"

namespace Mantid
{
namespace Algorithms
{

// Register the class into the algorithm factory
DECLARE_ALGORITHM(CalculateEfficiency)


using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;

/** Initialization method.
 *
 */
void CalculateEfficiency::init()
{
  declareProperty(
    new WorkspaceProperty<>("InputWorkspace","",Direction::Input),
                            "The workspace containing the flood data" );
  declareProperty(
    new WorkspaceProperty<>("OutputWorkspace","",Direction::Output),
    "The name of the workspace to be created as the output of the algorithm" );

  auto positiveDouble = boost::make_shared<BoundedValidator<double> >();
  positiveDouble->setLower(0);
  declareProperty("MinEfficiency", EMPTY_DBL(), positiveDouble,
      "Minimum efficiency for a pixel to be considered (default: no minimum).");
  declareProperty("MaxEfficiency", EMPTY_DBL(), positiveDouble,
      "Maximum efficiency for a pixel to be considered (default: no maximum).");

}

/** Executes the algorithm
 *
 */
void CalculateEfficiency::exec()
{

  // Minimum efficiency. Pixels with lower efficiency will be masked
  double min_eff = getProperty("MinEfficiency");
  // Maximum efficiency. Pixels with higher efficiency will be masked
  double max_eff = getProperty("MaxEfficiency");

  // Get the input workspace
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr rebinnedWS;// = inputWS;

  // Now create the output workspace
  MatrixWorkspace_sptr outputWS;// = getProperty("OutputWorkspace");

 // DataObjects::EventWorkspace_const_sptr inputEventWS = boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);

  // Sum up all the wavelength bins
  IAlgorithm_sptr childAlg = createChildAlgorithm("Integration", 0.0, 0.2);
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
  childAlg->executeAsChildAlg();
  rebinnedWS = childAlg->getProperty("OutputWorkspace");

  outputWS = WorkspaceFactory::Instance().create(rebinnedWS);
  WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, false);
  for (int i=0; i<(int)rebinnedWS->getNumberHistograms(); i++)
  {
    outputWS->dataX(i) = rebinnedWS->readX(i);
  }
  setProperty("OutputWorkspace",outputWS);

  double sum = 0.0;
  double err = 0.0;
  int npixels = 0;

  // Loop over spectra and sum all the counts to get normalization
  // Skip monitors and masked detectors
  sumUnmaskedDetectors(rebinnedWS, sum, err, npixels);

  // Normalize each detector pixel by the sum we just found to get the
  // relative efficiency. If the minimum and maximum efficiencies are
  // provided, the pixels falling outside this range will be marked
  // as 'masked' in both the input and output workspace.
  // We mask detectors in the input workspace so that we can resum the
  // counts to find a new normalization factor that takes into account
  // the newly masked detectors.
  normalizeDetectors(rebinnedWS, outputWS, sum, err, npixels, min_eff, max_eff);

  if ( !isEmpty(min_eff) || !isEmpty(max_eff) )
  {
    // Recompute the normalization, excluding the pixels that were outside
    // the acceptable efficiency range.
    sumUnmaskedDetectors(rebinnedWS, sum, err, npixels);

    // Now that we have a normalization factor that excludes bad pixels,
    // recompute the relative efficiency.
    // We pass EMPTY_DBL() to avoid masking pixels that might end up high or low
    // after the new normalization.
    normalizeDetectors(rebinnedWS, outputWS, sum, err, npixels, EMPTY_DBL(), EMPTY_DBL());
  }

  return;
}

/*
 *  Sum up all the unmasked detector pixels.
 *
 * @param rebinnedWS: workspace where all the wavelength bins have been grouped together
 * @param sum: sum of all the unmasked detector pixels (counts)
 * @param error: error on sum (counts)
 * @param nPixels: number of unmasked detector pixels that contributed to sum
 */
void CalculateEfficiency::sumUnmaskedDetectors(MatrixWorkspace_sptr rebinnedWS,
    double& sum, double& error, int& nPixels)
{
    // Number of spectra
    const int numberOfSpectra = static_cast<int>(rebinnedWS->getNumberHistograms());
    sum = 0.0;
    error = 0.0;
    nPixels = 0;

    for (int i = 0; i < numberOfSpectra; i++)
    {
      progress(0.2+0.2*i/numberOfSpectra, "Computing sensitivity");
      // Get the detector object for this spectrum
      IDetector_const_sptr det = rebinnedWS->getDetector(i);
      // If this detector is masked, skip to the next one
      if ( det->isMasked() ) continue;
      // If this detector is a monitor, skip to the next one
      if ( det->isMonitor() ) continue;

      // Retrieve the spectrum into a vector
      const MantidVec& YValues = rebinnedWS->readY(i);
      const MantidVec& YErrors = rebinnedWS->readE(i);

      sum += YValues[0];
      error += YErrors[0]*YErrors[0];
      nPixels++;
    }

    error = std::sqrt(error);
}

/*
 * Normalize each detector to produce the relative detector efficiency.
 * Pixels that fall outside those efficiency limits will be masked in both
 * the input and output workspaces.
 *
 * @param rebinnedWS: input workspace
 * @param outputWS: output workspace containing the relative efficiency
 * @param sum: sum of all the unmasked detector pixels (counts)
 * @param error: error on sum (counts)
 * @param nPixels: number of unmasked detector pixels that contributed to sum
 */

void CalculateEfficiency::normalizeDetectors(MatrixWorkspace_sptr rebinnedWS,
    MatrixWorkspace_sptr outputWS, double sum, double error, int nPixels,
    double min_eff, double max_eff)
{
    // Number of spectra
    const size_t numberOfSpectra = rebinnedWS->getNumberHistograms();

    // Empty vector to store the pixels that outside the acceptable efficiency range
    std::vector<size_t> dets_to_mask;

    for (size_t i = 0; i < numberOfSpectra; i++)
    {
      const double currProgress = 0.4+0.2*((double)i/(double)numberOfSpectra);
      progress(currProgress, "Computing sensitivity");
      // Get the detector object for this spectrum
      IDetector_const_sptr det = rebinnedWS->getDetector(i);
      // If this detector is masked, skip to the next one
      if ( det->isMasked() ) continue;

      // Retrieve the spectrum into a vector
      const MantidVec& YIn = rebinnedWS->readY(i);
      const MantidVec& EIn = rebinnedWS->readE(i);
      MantidVec& YOut = outputWS->dataY(i);
      MantidVec& EOut = outputWS->dataE(i);
      // If this detector is a monitor, skip to the next one
      if ( det->isMonitor() )
      {
        YOut[0] = 1.0;
        EOut[0] = 0.0;
        continue;
      }

      // Normalize counts to get relative efficiency
      YOut[0] = nPixels/sum * YIn[0];
      const double err_sum = YIn[0]/sum*error;
      EOut[0] = nPixels/std::abs(sum) * std::sqrt(EIn[0]*EIn[0] + err_sum*err_sum);

      // Mask this detector if the signal is outside the acceptable band
      if ( !isEmpty(min_eff) && YOut[0] < min_eff ) dets_to_mask.push_back(i);
      if ( !isEmpty(max_eff) && YOut[0] > max_eff ) dets_to_mask.push_back(i);

    }

    // If we identified pixels to be masked, mask them now
    if ( !dets_to_mask.empty() )
    {
      // Mask detectors that were found to be outside the acceptable efficiency band
      try
      {
        IAlgorithm_sptr mask = createChildAlgorithm("MaskDetectors", 0.8, 0.9);
        // First we mask detectors in the output workspace
        mask->setProperty<MatrixWorkspace_sptr>("Workspace", outputWS);
        mask->setProperty< std::vector<size_t> >("WorkspaceIndexList", dets_to_mask);
        mask->execute();

        mask = createChildAlgorithm("MaskDetectors", 0.9, 1.0);
        // Then we mask the same detectors in the input workspace
        mask->setProperty<MatrixWorkspace_sptr>("Workspace", rebinnedWS);
        mask->setProperty< std::vector<size_t> >("WorkspaceIndexList", dets_to_mask);
        mask->execute();
      } catch (std::invalid_argument& err)
      {
        std::stringstream e;
        e << "Invalid argument to MaskDetectors Child Algorithm: " << err.what();
        g_log.error(e.str());
      } catch (std::runtime_error& err)
      {
        std::stringstream e;
        e << "Unable to successfully run MaskDetectors Child Algorithm: " << err.what();
        g_log.error(e.str());
      }
    }
}

} // namespace Algorithms
} // namespace Mantid


syntax highlighted by Code2HTML, v. 0.9.1