MuonGroupAsymmetryCalculator.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 "MantidWorkflowAlgorithms/MuonGroupAsymmetryCalculator.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/VectorHelper.h"
 
using Mantid::API::AlgorithmManager;
using Mantid::API::IAlgorithm_sptr;
using Mantid::API::ITableWorkspace_sptr;
using Mantid::API::MatrixWorkspace_sptr;
using Mantid::API::Workspace_sptr;
using Mantid::API::WorkspaceFactory;
 
namespace Mantid::WorkflowAlgorithms {
 
//----------------------------------------------------------------------------------------------
MuonGroupAsymmetryCalculator::MuonGroupAsymmetryCalculator(const Mantid::API::WorkspaceGroup_sptr &inputWS,
                                                           const std::vector<int> &summedPeriods,
                                                           const std::vector<int> &subtractedPeriods,
                                                           const int groupIndex, const double start, const double end,
                                                           const std::string &wsName)
    : MuonGroupCalculator(inputWS, summedPeriods, subtractedPeriods, groupIndex) {
  MuonGroupCalculator::setStartEnd(start, end);
  MuonGroupCalculator::setWSName(wsName);
}
 
MatrixWorkspace_sptr MuonGroupAsymmetryCalculator::calculate() const {
  // The output workspace
 
  MatrixWorkspace_sptr tempWS;
  int numPeriods = m_inputWS->getNumberOfEntries();
  if (numPeriods > 1) {
    // Several period workspaces were supplied
 
    auto summedWS = sumPeriods(m_summedPeriods);
    auto subtractedWS = sumPeriods(m_subtractedPeriods);
 
    // Remove decay (summed periods ws)
    MatrixWorkspace_sptr asymSummedPeriods = estimateAsymmetry(summedWS, m_groupIndex);
 
    if (!m_subtractedPeriods.empty()) {
      // Remove decay (subtracted periods ws)
      MatrixWorkspace_sptr asymSubtractedPeriods = estimateAsymmetry(subtractedWS, m_groupIndex);
 
      // Now subtract
      tempWS = subtractWorkspaces(asymSummedPeriods, asymSubtractedPeriods);
    } else {
      tempWS = asymSummedPeriods;
    }
  } else {
    // Only one period was supplied
    tempWS = estimateAsymmetry(m_inputWS->getItem(0),
                               m_groupIndex); // change -1 to m_groupIndex and
                                              // follow through to store as a
                                              // table for later.
  }
 
  // Extract the requested spectrum
  MatrixWorkspace_sptr outWS = extractSpectrum(tempWS, m_groupIndex);
  return outWS;
}
 
MatrixWorkspace_sptr MuonGroupAsymmetryCalculator::removeExpDecay(const Workspace_sptr &inputWS,
                                                                  const int index) const {
  MatrixWorkspace_sptr outWS;
  // Remove decay
  if (inputWS) {
    auto asym = AlgorithmManager::Instance().create("RemoveExpDecay");
    asym->setChild(true);
    asym->setProperty("InputWorkspace", inputWS);
    if (index > 0) {
      // GroupIndex as vector
      // Necessary if we want RemoveExpDecay to fit only the requested
      // spectrum
      std::vector<int> spec(1, index);
      asym->setProperty("Spectra", spec);
    }
    asym->setProperty("OutputWorkspace", "__NotUsed__");
    asym->execute();
    outWS = asym->getProperty("OutputWorkspace");
  }
  return outWS;
}
MatrixWorkspace_sptr MuonGroupAsymmetryCalculator::estimateAsymmetry(const Workspace_sptr &inputWS,
                                                                     const int index) const {
  const std::string normTableName = "MuonAnalysisTFNormalizations";
  API::AnalysisDataServiceImpl &ads = API::AnalysisDataService::Instance();
  if (!ads.doesExist(normTableName)) {
    Mantid::API::ITableWorkspace_sptr table = Mantid::API::WorkspaceFactory::Instance().createTable();
    Mantid::API::AnalysisDataService::Instance().addOrReplace(normTableName, table);
    table->addColumn("double", "norm");
    table->addColumn("str", "name");
    table->addColumn("str", "method");
  }
 
  MatrixWorkspace_sptr outWS;
  // calculate asymmetry
  if (inputWS) {
    auto asym = AlgorithmManager::Instance().create("EstimateMuonAsymmetryFromCounts");
    asym->setChild(true);
 
    asym->setProperty("InputWorkspace", inputWS);
    asym->setProperty("WorkspaceName", m_wsName);
    if (index > -1) {
      std::vector<int> spec(1, index);
      asym->setProperty("Spectra", spec);
    }
    asym->setProperty("OutputWorkspace", "__NotUsed__");
    asym->setProperty("StartX", m_startX);
    asym->setProperty("EndX", m_endX);
    asym->setProperty("NormalizationIn", getStoredNorm());
    asym->setProperty("OutputUnNormData", true);
    asym->setProperty("OutputUnNormWorkspace", "tmp_unNorm");
    asym->setProperty("NormalizationTable", normTableName);
    asym->execute();
 
    API::MatrixWorkspace_sptr unnorm = asym->getProperty("OutputUnNormWorkspace");
    MatrixWorkspace_sptr singleWS = extractSpectrum(unnorm, index);
    ads.addOrReplace("tmp_unNorm", singleWS);
 
    outWS = asym->getProperty("OutputWorkspace");
  }
  return outWS;
}
/*
 * Reads in the stored normalization
 * if >0 then it is used instead of
 * estimating the norm.
 */
double getStoredNorm() {
  if (!API::AnalysisDataService::Instance().doesExist("__keepNorm__")) {
    return 0.0;
  } else {
    API::ITableWorkspace_sptr table =
        std::dynamic_pointer_cast<API::ITableWorkspace>(API::AnalysisDataService::Instance().retrieve("__keepNorm__"));
    auto colNorm = table->getColumn("norm");
    return (*colNorm)[0];
  }
}
 
} // namespace Mantid::WorkflowAlgorithms