ApplyDeadTimeCorr.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/ApplyDeadTimeCorr.h"
#include "MantidKernel/System.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/PropertyWithValue.h"
#include "boost/lexical_cast.hpp"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/TableRow.h"

#include <cmath>

using std::string;
using namespace Mantid::Kernel;
using namespace Mantid::API;

namespace Mantid {
namespace Algorithms {

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

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
   */
void ApplyDeadTimeCorr::init() {

  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>(
                      "InputWorkspace", "", Direction::Input),
                  "The name of the input workspace containing measured counts");

  declareProperty(new API::WorkspaceProperty<API::ITableWorkspace>(
                      "DeadTimeTable", "", Direction::Input),
                  "Name of the Dead Time Table");

  declareProperty(
      new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace", "",
                                                       Direction::Output),
      "The name of the output workspace containing corrected counts");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
   */
void ApplyDeadTimeCorr::exec() {
  // Get pointers to the workspace and dead time table
  API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
  API::ITableWorkspace_sptr deadTimeTable = getProperty("DeadTimeTable");

  if (!(deadTimeTable->rowCount() > inputWs->getNumberHistograms())) {
    // Get number of good frames from Run object. This also serves as
    // a test to see if valid input workspace has been provided

    const API::Run &run = inputWs->run();
    if (run.hasProperty("goodfrm")) {
      double numGoodFrames =
          boost::lexical_cast<double>(run.getProperty("goodfrm")->value());

      if (numGoodFrames == 0) {
        throw std::runtime_error(
            "Number of good frames in the workspace is zero");
      }

      // Duplicate the input workspace. Only need to change Y values based on
      // dead time corrections
      IAlgorithm_sptr duplicate = createChildAlgorithm("CloneWorkspace");
      duplicate->initialize();
      duplicate->setProperty<Workspace_sptr>(
          "InputWorkspace", boost::dynamic_pointer_cast<Workspace>(inputWs));
      duplicate->execute();
      Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
      MatrixWorkspace_sptr outputWs =
          boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

      // Presumed to be the same for all data
      double timeBinWidth(inputWs->dataX(0)[1] - inputWs->dataX(0)[0]);

      if (timeBinWidth != 0) {
        try {
          // Apply Dead Time
          for (size_t i = 0; i < deadTimeTable->rowCount(); ++i) {
            API::TableRow deadTimeRow = deadTimeTable->getRow(i);
            size_t index =
                static_cast<size_t>(inputWs->getIndexFromSpectrumNumber(
                    static_cast<int>(deadTimeRow.Int(0))));

            for (size_t j = 0; j < inputWs->blocksize(); ++j) {
              double temp(
                  1 -
                  inputWs->dataY(index)[j] *
                      (deadTimeRow.Double(1) / (timeBinWidth * numGoodFrames)));
              if (temp != 0) {
                outputWs->dataY(index)[j] = inputWs->dataY(index)[j] / temp;
              } else {
                g_log.error() << "1 - MeasuredCount * (Deadtime/TimeBin width "
                                 "is currently (" << temp
                              << "). Can't divide by this amount." << std::endl;

                throw std::invalid_argument("Can't divide by 0");
              }
            }
          }

          setProperty("OutputWorkspace", outputWs);
        } catch (std::runtime_error &) {
          throw std::invalid_argument("Invalid argument for algorithm.");
        }
      } else {
        g_log.error() << "The time bin width is currently (" << timeBinWidth
                      << "). Can't divide by this amount." << std::endl;

        throw std::invalid_argument("Can't divide by 0");
      }
    } else {
      g_log.error() << "To calculate Muon deadtime requires that goodfrm "
                       "(number of good frames) "
                    << "is stored in InputWorkspace Run object\n";
    }
  } else {
    g_log.error()
        << "Row count(" << deadTimeTable->rowCount()
        << ") of Dead time table is bigger than the Number of Histograms("
        << inputWs->getNumberHistograms() << ")." << std::endl;

    throw std::invalid_argument(
        "Row count was bigger than the Number of Histograms.");
  }
}

} // namespace Mantid
} // namespace Algorithms