EditInstrumentGeometry.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/EditInstrumentGeometry.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/MandatoryValidator.h"
 
#include <sstream>
 
using namespace Mantid::Kernel;
using namespace Mantid::API;
 
using namespace std;
 
namespace Mantid::Algorithms {
 
DECLARE_ALGORITHM(EditInstrumentGeometry)
 
const std::string EditInstrumentGeometry::name() const { return "EditInstrumentGeometry"; }
 
const std::string EditInstrumentGeometry::category() const {
  return "Diffraction\\DataHandling;DataHandling\\Instrument";
}
 
int EditInstrumentGeometry::version() const { return 1; }
 
//----------------------------------------------------------------------------------------------
 
//----------------------------------------------------------------------------------------------
void EditInstrumentGeometry::init() {
  // Input workspace
  declareProperty(std::make_unique<WorkspaceProperty<>>("Workspace", "", Direction::InOut),
                  "Workspace to edit the detector information");
 
  // L1
  declareProperty("PrimaryFlightPath", EMPTY_DBL(), "Primary flight path L1 of the powder diffractometer. ");
 
  // Spectrum Number for the spectrum to have instrument geometry edited
  declareProperty(std::make_unique<ArrayProperty<int32_t>>("SpectrumIDs"),
                  "Spectrum Numbers (note that it is not detector ID or workspace "
                  "indices). The list must be either empty or have a size "
                  "equal to input workspace's histogram number. ");
 
  auto required = std::make_shared<MandatoryValidator<std::vector<double>>>();
 
  // Vector for L2
  declareProperty(std::make_unique<ArrayProperty<double>>("L2", required),
                  "Secondary flight (L2) paths for each detector.  Number of L2 "
                  "given must be same as number of histogram.");
 
  // Vector for 2Theta
  declareProperty(std::make_unique<ArrayProperty<double>>("Polar", required),
                  "Polar angles (two thetas) for detectors. Number of 2theta "
                  "given must be same as number of histogram.");
 
  // Vector for Azimuthal angle
  declareProperty(std::make_unique<ArrayProperty<double>>("Azimuthal"),
                  "Azimuthal angles (out-of-plane) for detectors. "
                  "Number of azimuthal angles given must be same as number of histogram.");
 
  // Detector IDs
  declareProperty(std::make_unique<ArrayProperty<int>>("DetectorIDs"),
                  "User specified detector IDs of the spectra. "
                  "Number of specified detector IDs must be either zero or "
                  "number of histogram");
 
  // Instrument Name
  declareProperty("InstrumentName", "",
                  "Name of the newly built instrument.  If left empty, "
                  "the original instrument will be used. ");
}
 
template <typename NumT> std::string checkValues(const std::vector<NumT> &thingy, const size_t numHist) {
  if ((!thingy.empty()) && thingy.size() != numHist) {
    stringstream msg;
    msg << "Must equal number of spectra or be empty (" << numHist << " != " << thingy.size() << ")";
    return msg.str();
  } else {
    return "";
  }
}
 
std::map<std::string, std::string> EditInstrumentGeometry::validateInputs() {
  std::map<std::string, std::string> result;
 
  // everything depends on being parallel to the workspace # histo
  size_t numHist(0);
  bool hasWorkspacePtr(false);
  {
    MatrixWorkspace_const_sptr workspace = getProperty("Workspace");
    // this is to guard against WorkspaceGroups
    // fallthrough is to skip workspace check and make sure
    if (bool(workspace)) {
      hasWorkspacePtr = true;
      numHist = workspace->getNumberHistograms();
    }
  }
 
  std::string error;
 
  const std::vector<int32_t> specids = this->getProperty("SpectrumIDs");
  if (!hasWorkspacePtr) {
    // use the number of spectra for the number of histograms
    numHist = specids.size();
    // give up if it is empty
    if (numHist == 0) {
      return result;
    }
  }
  error = checkValues(specids, numHist);
  if (!error.empty())
    result["SpectrumIDs"] = error;
 
  const std::vector<double> l2 = this->getProperty("L2");
  error = checkValues(l2, numHist);
  if (!error.empty())
    result["L2"] = error;
 
  const std::vector<double> tth = this->getProperty("Polar");
  error = checkValues(tth, numHist);
  if (!error.empty())
    result["Polar"] = error;
 
  const std::vector<double> phi = this->getProperty("Azimuthal");
  error = checkValues(phi, numHist);
  if (!error.empty())
    result["Azimuthal"] = error;
 
  const vector<int> detids = getProperty("DetectorIDs");
  error = checkValues(detids, numHist);
  if (!error.empty())
    result["DetectorIDs"] = error;
 
  // TODO verify that SpectrumIDs, L2, Polar, Azimuthal, and DetectorIDs are
  // parallel or not specified
  return result;
}
 
//----------------------------------------------------------------------------------------------
void EditInstrumentGeometry::exec() {
  // Lots of things have to do with the input workspace
  MatrixWorkspace_sptr workspace = getProperty("Workspace");
  Geometry::Instrument_const_sptr originstrument = workspace->getInstrument();
 
  // Get and check the primary flight path
  double l1 = this->getProperty("PrimaryFlightPath");
  if (isEmpty(l1)) {
    // Use the original L1
    if (!originstrument) {
      std::string errmsg("It is not supported that L1 is not given, ",
                         "while there is no instrument associated to input workspace.");
      g_log.error(errmsg);
      throw std::runtime_error(errmsg);
    }
    Geometry::IComponent_const_sptr source = originstrument->getSource();
    Geometry::IComponent_const_sptr sample = originstrument->getSample();
    l1 = source->getDistance(*sample);
    g_log.information() << "Retrieve L1 from input data workspace. \n";
  }
  g_log.information() << "Using L1 = " << l1 << "\n";
 
  // Get spectra number in case they are in a funny order
  std::vector<int32_t> specids = this->getProperty("SpectrumIDs");
  if (specids.empty()) // they are using the order of the input workspace
  {
    size_t numHist = workspace->getNumberHistograms();
    for (size_t i = 0; i < numHist; ++i) {
      specids.emplace_back(workspace->getSpectrum(i).getSpectrumNo());
      g_log.information() << "Add spectrum " << workspace->getSpectrum(i).getSpectrumNo() << ".\n";
    }
  }
 
  // Get the detector ids - empsy means ignore it
  const vector<int> vec_detids = getProperty("DetectorIDs");
  const bool renameDetID(!vec_detids.empty());
 
  // Get individual detector geometries ordered by input spectrum Numbers
  const std::vector<double> l2s = this->getProperty("L2");
  const std::vector<double> tths = this->getProperty("Polar");
  std::vector<double> phis = this->getProperty("Azimuthal");
 
  // empty list of L2 and 2-theta is not allowed
  if (l2s.empty()) {
    throw std::runtime_error("User must specify L2 for all spectra. ");
  }
  if (tths.empty()) {
    throw std::runtime_error("User must specify 2theta for all spectra.");
  }
 
  // empty list of phi means that they are all zero
  if (phis.empty()) {
    phis.assign(l2s.size(), 0.);
  }
 
  // Validate
  for (size_t ib = 0; ib < l2s.size(); ib++) {
    g_log.information() << "Detector " << specids[ib] << "  L2 = " << l2s[ib] << "  2Theta = " << tths[ib] << '\n';
    if (specids[ib] < 0) {
      // Invalid spectrum Number : less than 0.
      stringstream errmsgss;
      errmsgss << "Detector ID = " << specids[ib] << " cannot be less than 0.";
      throw std::invalid_argument(errmsgss.str());
    }
    if (l2s[ib] <= 0.0) {
      throw std::invalid_argument("L2 cannot be less or equal to 0");
    }
  }
 
  // Keep original instrument and set the new instrument, if necessary
  const auto spec2indexmap = workspace->getSpectrumToWorkspaceIndexMap();
 
  // ??? Condition: spectrum has 1 and only 1 detector
  size_t nspec = workspace->getNumberHistograms();
 
  // Initialize another set of L2/2-theta/Phi/DetectorIDs vector ordered by
  // workspace index
  std::vector<double> storL2s(nspec, 0.);
  std::vector<double> stor2Thetas(nspec, 0.);
  std::vector<double> storPhis(nspec, 0.);
  vector<int> storDetIDs(nspec, 0);
 
  // Map the properties from spectrum Number to workspace index
  for (size_t i = 0; i < specids.size(); i++) {
    // Find spectrum's workspace index
    auto it = spec2indexmap.find(specids[i]);
    if (it == spec2indexmap.end()) {
      stringstream errss;
      errss << "Spectrum Number " << specids[i] << " is not found. "
            << "Instrument won't be edited for this spectrum. \n";
      g_log.error(errss.str());
      throw std::runtime_error(errss.str());
    }
 
    // Store and set value
    size_t workspaceindex = it->second;
 
    storL2s[workspaceindex] = l2s[i];
    stor2Thetas[workspaceindex] = tths[i];
    storPhis[workspaceindex] = phis[i];
    if (renameDetID)
      storDetIDs[workspaceindex] = vec_detids[i];
 
    g_log.debug() << "workspace index = " << workspaceindex << " is for Spectrum " << specids[i] << '\n';
  }
 
  // Generate a new instrument
  // Name of the new instrument
  std::string name = std::string(getProperty("InstrumentName"));
  if (name.empty()) {
    // Use the original L1
    if (!originstrument) {
      std::string errmsg("It is not supported that InstrumentName is not given, ",
                         "while there is no instrument associated to input workspace.");
      g_log.error(errmsg);
      throw std::runtime_error(errmsg);
    }
    name = originstrument->getName();
  }
 
  // Create a new instrument from scratch any way.
  auto instrument = std::make_shared<Geometry::Instrument>(name);
  if (!bool(instrument)) {
    stringstream errss;
    errss << "Trying to use a Parametrized Instrument as an Instrument.";
    g_log.error(errss.str());
    throw std::runtime_error(errss.str());
  }
 
  // Set up source and sample information
  Geometry::Component *samplepos = new Geometry::Component("Sample", instrument.get());
  instrument->add(samplepos);
  instrument->markAsSamplePos(samplepos);
  samplepos->setPos(0.0, 0.0, 0.0);
 
  Geometry::ObjComponent *source = new Geometry::ObjComponent("Source", instrument.get());
  instrument->add(source);
  instrument->markAsSource(source);
  source->setPos(0.0, 0.0, -1.0 * l1);
 
  // Add/copy detector information
  for (size_t i = 0; i < workspace->getNumberHistograms(); i++) {
    // Create a new detector.
    //    (Instrument will take ownership of pointer so no need to delete.)
    detid_t newdetid;
    if (renameDetID)
      newdetid = storDetIDs[i];
    else
      newdetid = detid_t(i) + 100;
    Geometry::Detector *detector = new Geometry::Detector("det", newdetid, samplepos);
 
    // Set up new detector parameters related to new instrument
    double l2 = storL2s[i];
    double tth = stor2Thetas[i];
    double phi = storPhis[i];
 
    Kernel::V3D pos;
    pos.spherical(l2, tth, phi);
    detector->setPos(pos);
 
    // Add new detector to spectrum and instrument
    auto &spectrum = workspace->getSpectrum(i);
    // Good and do some debug output
    g_log.debug() << "Orignal spectrum " << spectrum.getSpectrumNo() << "has " << spectrum.getDetectorIDs().size()
                  << " detectors. \n";
 
    spectrum.clearDetectorIDs();
    spectrum.addDetectorID(newdetid);
    instrument->add(detector);
    instrument->markAsDetector(detector);
 
  } // ENDFOR workspace index
 
  // Add the new instrument
  workspace->setInstrument(instrument);
}
 
} // namespace Mantid::Algorithms