RemovePromptPulse.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/RemovePromptPulse.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/TimeSeriesProperty.h"
 
using std::string;
namespace Mantid::Algorithms {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(RemovePromptPulse)
 
using namespace Mantid::Kernel;
using namespace Mantid::API;
 
//----------------------------------------------------------------------------------------------
 
const string RemovePromptPulse::name() const { return "RemovePromptPulse"; }
 
int RemovePromptPulse::version() const { return 1; }
 
const string RemovePromptPulse::category() const { return "CorrectionFunctions\\BackgroundCorrections"; }
 
//----------------------------------------------------------------------------------------------
void RemovePromptPulse::init() {
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
                                                        std::make_shared<WorkspaceUnitValidator>("TOF")),
                  "An input workspace.");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "An output workspace.");
 
  auto validator = std::make_shared<BoundedValidator<double>>();
  validator->setLower(0.0);
  declareProperty("Width", Mantid::EMPTY_DBL(), validator,
                  "The width of the time of flight (in microseconds) to remove "
                  "from the data.");
  declareProperty("Frequency", Mantid::EMPTY_DBL(), validator,
                  "The frequency of the source (in Hz) used to calculate the "
                  "minimum time of flight to filter.");
}
 
//----------------------------------------------------------------------------------------------
namespace { // anonymous namespace begin
double getMedian(const API::Run &run, const std::string &name) {
 
  if (!run.hasProperty(name)) {
    return Mantid::EMPTY_DBL();
  }
  auto *log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(run.getLogData(name));
  if (!log)
    return Mantid::EMPTY_DBL();
 
  Kernel::TimeSeriesPropertyStatistics stats = log->getStatistics();
  return stats.median;
}
 
void getTofRange(const MatrixWorkspace_const_sptr &wksp, double &tmin, double &tmax) {
  DataObjects::EventWorkspace_const_sptr eventWksp = std::dynamic_pointer_cast<const DataObjects::EventWorkspace>(wksp);
  if (eventWksp == nullptr) {
    wksp->getXMinMax(tmin, tmax);
  } else {
    eventWksp->getEventXMinMax(tmin, tmax);
  }
}
} // namespace
 
void RemovePromptPulse::exec() {
  // verify there is a width parameter specified
  double width = this->getProperty("Width");
  if (this->isEmpty(width)) {
    throw std::runtime_error("Failed to specify \'Width\' parameter");
  }
 
  // need the input workspace in general
  API::MatrixWorkspace_const_sptr inputWS = this->getProperty("InputWorkspace");
 
  // get the frequency
  double frequency = this->getProperty("Frequency");
  if (this->isEmpty(frequency)) // it wasn't specified so try divination
  {
    frequency = this->getFrequency(inputWS->run());
    if (this->isEmpty(frequency)) {
      throw std::runtime_error("Failed to determine the frequency");
    }
  }
  g_log.information() << "Using frequency of " << frequency << "Hz\n";
  double period = 1000000. / frequency; // period in microseconds
 
  // determine the overall tof window for the data
  double tmin;
  double tmax;
  getTofRange(inputWS, tmin, tmax);
  g_log.information() << "Data tmin=" << tmin << ", tmax=" << tmax << ", period=" << period << " microseconds\n";
 
  // calculate the times for the prompt pulse
  std::vector<double> pulseTimes = this->calculatePulseTimes(tmin, tmax, period);
  if (pulseTimes.empty()) {
    g_log.notice() << "Not applying filter since prompt pulse is not in data "
                      "range (period = "
                   << period << ")\n";
    setProperty("OutputWorkspace", std::const_pointer_cast<MatrixWorkspace>(inputWS));
    return;
  }
  g_log.information() << "Calculated prompt pulses at ";
  for (double pulseTime : pulseTimes)
    g_log.information() << pulseTime << " ";
  g_log.information() << " microseconds\n";
 
  MatrixWorkspace_sptr outputWS;
  for (double &pulseTime : pulseTimes) {
    double right = pulseTime + width;
 
    g_log.notice() << "Filtering tmin=" << pulseTime << ", tmax=" << right << " microseconds\n";
 
    // run maskbins to do the work on the first prompt pulse
    auto algo = createChildAlgorithm("MaskBins");
    if (outputWS) {
      algo->setProperty<MatrixWorkspace_sptr>("InputWorkspace", std::const_pointer_cast<MatrixWorkspace>(outputWS));
    } else { // should only be first time
      algo->setProperty<MatrixWorkspace_sptr>("InputWorkspace", std::const_pointer_cast<MatrixWorkspace>(inputWS));
      outputWS = this->getProperty("OutputWorkspace");
    }
    // always write to correct output workspace
    algo->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputWS);
    algo->setProperty<double>("XMin", pulseTime);
    algo->setProperty<double>("XMax", right);
    algo->executeAsChildAlg();
 
    // copy over the output workspace
    outputWS = algo->getProperty("OutputWorkspace");
  }
 
  setProperty("OutputWorkspace", outputWS);
}
 
double RemovePromptPulse::getFrequency(const API::Run &run) {
  double candidate;
 
  candidate = getMedian(run, "Frequency");
  if (!(this->isEmpty(candidate)))
    return candidate;
 
  candidate = getMedian(run, "frequency");
  if (!(this->isEmpty(candidate)))
    return candidate;
 
  candidate = getMedian(run, "FREQUENCY");
  if (!(this->isEmpty(candidate)))
    return candidate;
 
  // give up
  return Mantid::EMPTY_DBL();
}
 
std::vector<double> RemovePromptPulse::calculatePulseTimes(const double tmin, const double tmax, const double period) {
  std::vector<double> times;
  double time = 0.;
 
  // find when the first prompt pulse would be
  while (time < tmin)
    time += period;
 
  // calculate all times possible
  while (time < tmax) {
    times.emplace_back(time);
    time += period;
  }
 
  return times;
}
} // namespace Mantid::Algorithms