IntegrateByComponent.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/IntegrateByComponent.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidTypes/SpectrumDefinition.h"
 
#include <gsl/gsl_statistics.h>
#include <unordered_map>
 
namespace Mantid::Algorithms {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IntegrateByComponent)
 
using namespace Mantid::API;
using namespace Mantid::Kernel;
 
//----------------------------------------------------------------------------------------------
const std::string IntegrateByComponent::name() const { return "IntegrateByComponent"; }
 
int IntegrateByComponent::version() const { return 1; }
 
const std::string IntegrateByComponent::category() const { return "Utility\\Workspaces"; }
 
//----------------------------------------------------------------------------------------------
 
//----------------------------------------------------------------------------------------------
void IntegrateByComponent::init() {
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
                                                        std::make_shared<HistogramValidator>()),
                  "The input workspace.");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "The output workspace.");
  auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
  mustBePosInt->setLower(0);
  declareProperty("LevelsUp", 0, mustBePosInt,
                  "Levels above pixel that will be used to compute the average.\n"
                  "If no level is specified, the median is over the whole instrument.\n If "
                  "0, it will just return the integrated values in each pixel");
}
 
//----------------------------------------------------------------------------------------------
void IntegrateByComponent::exec() {
  MatrixWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
  int parents = getProperty("LevelsUp");
  // Make sure it's integrated
  auto childAlg = createChildAlgorithm("Integration", 0, 0.2);
  childAlg->setProperty("InputWorkspace", inputWS);
  childAlg->setProperty("StartWorkspaceIndex", 0);
  childAlg->setProperty("EndWorkspaceIndex", EMPTY_INT());
  childAlg->setProperty("RangeLower", inputWS->getXMin());
  childAlg->setProperty("RangeUpper", inputWS->getXMax());
  childAlg->setPropertyValue("IncludePartialBins", "1");
  childAlg->executeAsChildAlg();
  MatrixWorkspace_sptr integratedWS = childAlg->getProperty("OutputWorkspace");
 
  if (parents > 0) {
    std::vector<std::vector<size_t>> specmap = makeMap(integratedWS, parents);
    API::Progress prog(this, 0.3, 1.0, specmap.size());
    // calculate averages
    const auto &spectrumInfo = integratedWS->spectrumInfo();
    for (auto hists : specmap) {
      prog.report();
      std::vector<double> averageYInput, averageEInput;
 
      PARALLEL_FOR_IF(Kernel::threadSafe(*integratedWS))
      for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
        PARALLEL_START_INTERRUPT_REGION
 
        if (spectrumInfo.isMonitor(hists[i]))
          continue;
        if (spectrumInfo.isMasked(hists[i]))
          continue;
 
        const double yValue = integratedWS->y(hists[i])[0];
        const double eValue = integratedWS->e(hists[i])[0];
 
        if (!std::isfinite(yValue) || !std::isfinite(eValue)) // NaNs/Infs
          continue;
 
        // Now we have a good value
        PARALLEL_CRITICAL(IntegrateByComponent_good) {
          averageYInput.emplace_back(yValue);
          averageEInput.emplace_back(eValue * eValue);
        }
 
        PARALLEL_END_INTERRUPT_REGION
      }
      PARALLEL_CHECK_INTERRUPT_REGION
 
      double averageY, averageE;
      if (averageYInput.empty()) {
        g_log.information("some group has no valid histograms. Will use 0 for average.");
        averageY = 0.;
        averageE = 0.;
      } else {
        averageY = gsl_stats_mean(&averageYInput[0], 1, averageEInput.size());
        averageE = std::sqrt(gsl_stats_mean(&averageEInput[0], 1, averageYInput.size()));
      }
 
      PARALLEL_FOR_IF(Kernel::threadSafe(*integratedWS))
      for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
        PARALLEL_START_INTERRUPT_REGION
        if (spectrumInfo.isMonitor(hists[i]))
          continue;
        if (spectrumInfo.isMasked(hists[i]))
          continue;
 
        const double yValue = integratedWS->y(hists[i])[0];
        const double eValue = integratedWS->e(hists[i])[0];
        if (!std::isfinite(yValue) || !std::isfinite(eValue)) // NaNs/Infs
          continue;
 
        // Now we have a good value
        PARALLEL_CRITICAL(IntegrateByComponent_setaverage) {
          integratedWS->mutableY(hists[i])[0] = averageY;
          integratedWS->mutableE(hists[i])[0] = averageE;
        }
 
        PARALLEL_END_INTERRUPT_REGION
      }
      PARALLEL_CHECK_INTERRUPT_REGION
    }
  }
  // Assign it to the output workspace property
  setProperty("OutputWorkspace", integratedWS);
}
 
std::vector<std::vector<size_t>> IntegrateByComponent::makeInstrumentMap(const API::MatrixWorkspace_sptr &countsWS) {
  std::vector<std::vector<size_t>> mymap;
  std::vector<size_t> single;
 
  for (size_t i = 0; i < countsWS->getNumberHistograms(); i++) {
    single.emplace_back(i);
  }
  mymap.emplace_back(single);
  return mymap;
}
 
std::vector<std::vector<size_t>> IntegrateByComponent::makeMap(const API::MatrixWorkspace_sptr &countsWS, int parents) {
  // keyed by the component index of the ancestor
  std::unordered_multimap<size_t, size_t> mymap;
 
  if (parents == 0) // this should not happen in this file, but if one reuses
                    // the function and parents==0, the program has a sudden end
                    // without this check.
  {
    return makeInstrumentMap(countsWS);
  }
 
  const auto spectrumInfo = countsWS->spectrumInfo();
  const auto &componentInfo = countsWS->componentInfo();
  for (size_t i = 0; i < countsWS->getNumberHistograms(); i++) {
    if (!spectrumInfo.hasDetectors(i)) {
      g_log.debug("Spectrum has no detector, skipping");
      continue;
    }
 
    // walk up the component tree to the requested ancestor
    size_t ancestor = spectrumInfo.firstDetectorIndex(i);
    for (int level = 0; level < parents; ++level) {
      if (!componentInfo.hasParent(ancestor)) {
        g_log.warning("Too many levels up. Will ignore LevelsUp");
        return makeInstrumentMap(countsWS);
      }
      ancestor = componentInfo.parent(ancestor);
    }
    mymap.emplace(ancestor, i);
  }
 
  std::vector<std::vector<size_t>> speclist;
  std::vector<size_t> speclistsingle;
 
  auto s_it = mymap.begin();
  for (auto m_it = mymap.begin(); m_it != mymap.end(); m_it = s_it) {
    const size_t theKey = (*m_it).first;
    auto keyRange = mymap.equal_range(theKey);
 
    // Iterate over all map elements with key == theKey
    speclistsingle.clear();
    for (s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
      speclistsingle.emplace_back((*s_it).second);
    }
    speclist.emplace_back(speclistsingle);
  }
 
  return speclist;
}
 
} // namespace Mantid::Algorithms