SaveIsawPeaks.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 "MantidCrystal/SaveIsawPeaks.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidDataObjects/Peak.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/FacilityInfo.h"
#include "MantidKernel/InstrumentInfo.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Utils.h"
#include <boost/algorithm/string/trim.hpp>
#include <filesystem>
#include <fstream>
 
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
using namespace Mantid::Kernel;
using namespace Mantid::API;
 
namespace Mantid::Crystal {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SaveIsawPeaks)
 
 
void SaveIsawPeaks::init() {
  declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("InputWorkspace", "", Direction::Input,
                                                                      std::make_shared<InstrumentValidator>()),
                  "An input PeaksWorkspace with an instrument.");
 
  declareProperty("AppendFile", false,
                  "Append to file if true.\n"
                  "If false, new file (default).");
 
  const std::vector<std::string> exts{".peaks", ".integrate"};
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Save, exts),
                  "Path to an ISAW-style peaks or integrate file to save.");
 
  declareProperty(std::make_unique<WorkspaceProperty<Workspace2D>>("ProfileWorkspace", "", Direction::Input,
                                                                   PropertyMode::Optional),
                  "An optional Workspace2D of profiles from integrating cylinder.");
 
  declareProperty("RenumberPeaks", false,
                  "If true, sequential peak numbers\n"
                  "If false, keep original numbering (default).");
}
 
void SaveIsawPeaks::exec() {
  // Section header
  std::string header = "2   SEQN    H    K    L     COL      ROW     CHAN      "
                       "  L2   2_THETA        AZ         WL         D      IPK "
                       "      INTI    SIGI  RFLG";
 
  const std::string filename = getPropertyValue("Filename");
  m_ws = getProperty("InputWorkspace");
  const auto &peaks = m_ws->getPeaks();
  inst = m_ws->getInstrument();
  if (!inst)
    throw std::runtime_error("No instrument in the Workspace. Cannot save DetCal file.");
  const auto &detectorInfo = m_ws->detectorInfo();
 
  // We must sort the peaks first by run, then bank #, and save the list of
  // workspace indices of it
  using bankMap_t = std::map<int, std::vector<size_t>>;
  using runMap_t = std::map<int, bankMap_t>;
  std::set<int, std::less<int>> uniqueBanks;
  // We cannot assume the peaks have bank type detector modules, so we have a
  // string to check this
  std::string bankPart = "bank";
  if (inst->getName() == "WISH")
    bankPart = "WISHpanel";
 
  // Get all children
  const auto &componentInfo = m_ws->componentInfo();
  // iterate over the top level components, which contain the banks
  size_t const rootIndex = componentInfo.root();
  auto const topChildren = componentInfo.children(rootIndex);
  for (size_t const i : topChildren) {
    size_t bankParent = componentInfo.findBankParent(i, bankPart);
    auto const children = componentInfo.children(bankParent);
    for (size_t const child : children) {
      std::string bankName = componentInfo.name(child);
      boost::trim(bankName);
      boost::erase_all(bankName, bankPart);
      int bank = 0;
      Strings::convert(bankName, bank);
      if (bank == 0)
        continue;
      if (bankMasked(child, detectorInfo))
        continue;
      // Track unique bank numbers
      uniqueBanks.insert(bank);
    }
  }
  runMap_t runMap;
  for (size_t i = 0; i < peaks.size(); ++i) {
    const Peak &p = peaks[i];
    if (p.getIntMNP() != V3D(0, 0, 0))
      m_isModulatedStructure = true;
    int run = p.getRunNumber();
    int bank = 0;
    std::string bankName = p.getBankName();
    if (bankName.size() <= 4) {
      g_log.information() << "Could not interpret bank number of peak " << i << "(" << bankName << ")\n";
      continue;
    }
    // Save the "bank" part once to check whether it really is a bank
    if (bankPart == "?")
      bankPart = bankName.substr(0, 4);
    // Take out the "bank" part of the bank name and convert to an int
    if (bankPart == "bank")
      bankName = bankName.substr(4, bankName.size() - 4);
    else if (bankPart == "WISHpanel")
      bankName = bankName.substr(9, bankName.size() - 9);
    Strings::convert(bankName, bank);
 
    // Save in the map
    runMap[run][bank].emplace_back(i);
  }
  if (m_isModulatedStructure)
    header = "2   SEQN    H    K    L    M    N    P     COL      ROW     CHAN      "
             "  L2   2_THETA        AZ         WL         D      IPK "
             "      INTI    SIGI  RFLG";
 
  if (!inst)
    throw std::runtime_error("No instrument in PeaksWorkspace. Cannot save peaks file.");
 
  if (bankPart != "bank" && bankPart != "WISHpanel" && bankPart != "?") {
    std::ostringstream mess;
    mess << "Detector module of type " << bankPart << " not supported in ISAWPeaks. Cannot save peaks file";
    throw std::runtime_error(mess.str());
  }
 
  double l1;
  V3D beamline;
  double beamline_norm;
  V3D samplePos;
  inst->getInstrumentParameters(l1, beamline, beamline_norm, samplePos);
 
  std::ofstream out;
  bool append = getProperty("AppendFile");
  const bool renumber = getProperty("RenumberPeaks");
 
  // do not append if file does not exist
  if (!std::filesystem::exists(filename))
    append = false;
 
  int appendPeakNumb = 0;
  if (append) {
    std::ifstream infile(filename.c_str());
    std::string line;
    while (!infile.eof()) // To get you all the lines.
    {
      getline(infile, line); // Saves the line in STRING.
      if (infile.eof())
        break;
      std::stringstream ss(line);
      double three;
      ss >> three;
      if (three == 3) {
        int peakNumber;
        ss >> peakNumber;
        appendPeakNumb = std::max(peakNumber, appendPeakNumb);
      }
    }
 
    infile.close();
    out.open(filename.c_str(), std::ios::app);
    appendPeakNumb = appendPeakNumb + 1;
  } else {
    out.open(filename.c_str());
 
    const auto instrumentName = inst->getName();
    std::string facilityName;
    try {
      facilityName = ConfigService::Instance().getInstrument(instrumentName).facility().name();
    } catch (Exception::NotFoundError &) {
      g_log.warning() << "Instrument " << instrumentName
                      << " not found at any defined facility. Setting facility "
                         "name to Unknown\n";
      facilityName = "Unknown";
    }
    out << "Version: 2.0  Facility: " << facilityName;
    out << "  Instrument: " << instrumentName << "  Date: ";
 
    // TODO: The experiment date might be more useful than the instrument date.
    // For now, this allows the proper instrument to be loaded back after
    // saving.
    Types::Core::DateAndTime expDate = inst->getValidFromDate() + 1.0;
    out << expDate.toISO8601String();
    if (m_isModulatedStructure)
      out << " MOD";
    out << '\n';
 
    out << "6         L1    T0_SHIFT\n";
    out << "7 " << std::setw(10);
    out << std::setprecision(4) << std::fixed << (l1 * 100);
    out << std::setw(12) << std::setprecision(3) << std::fixed;
    // Time offset from property
    const API::Run &run = m_ws->run();
    double T0 = 0.0;
    if (run.hasProperty("T0")) {
      T0 = run.getPropertyValueAsType<double>("T0");
      if (T0 != 0) {
        g_log.notice() << "T0 = " << T0 << '\n';
      }
    }
    out << T0 << '\n';
 
    //  Save .detcal info
    out << "4 DETNUM  NROWS  NCOLS   WIDTH   HEIGHT   DEPTH   DETD   CenterX "
           "  CenterY   CenterZ    BaseX    BaseY    BaseZ      UpX      UpY "
           "     UpZ\n";
    // Here would save each detector...
    for (const auto bank : uniqueBanks) {
      // Build up the bank name
      std::ostringstream mess;
      if (bankPart == "bank")
        mess << "bank" << bank;
      else if (bankPart == "WISHpanel") {
        mess << "WISHpanel" << std::setfill('0') << std::setw(2) << bank;
      }
 
      std::string bankName = mess.str();
      // Retrieve it
      std::shared_ptr<const IComponent> det = inst->getComponentByName(bankName);
      if (inst->getName() == "CORELLI") // for Corelli with sixteenpack under bank
      {
        const size_t bankIndex = componentInfo.indexOfAny(bankName);
        const auto children = componentInfo.children(bankIndex);
        if (!children.empty()) {
          det.reset(componentInfo.componentID(children[0]), NoDeleting());
        }
      }
      if (det) {
        // Center of the detector.
        // NOTE: componentID() above returns the *base* component, whose getPos() ignores the ParameterMap and so
        //       reports uncalibrated IDF geometry. ComponentInfo is the parameterized view, so query it instead.
        V3D center = componentInfo.position(componentInfo.indexOf(det->getComponentID()));
 
        // Distance to center of detector
        double detd = (center - inst->getSample()->getPos()).norm();
        int NCOLS, NROWS;
        double xsize, ysize;
        sizeBanks(bankName, NCOLS, NROWS, xsize, ysize);
        // Base unit vector (along the horizontal, X axis)
        int midX = NCOLS / 2;
        int midY = NROWS / 2;
        V3D base = findPixelPos(bankName, midX + 1, midY) - findPixelPos(bankName, midX, midY);
        base.normalize();
 
        // Up unit vector (along the vertical, Y axis)
        V3D up = findPixelPos(bankName, midX, midY + 1) - findPixelPos(bankName, midX, midY);
        up.normalize();
 
        // Write the line
        out << "5 " << std::setw(6) << std::right << bank << " " << std::setw(6) << std::right << NROWS << " "
            << std::setw(6) << std::right << NCOLS << " " << std::setw(7) << std::right << std::fixed
            << std::setprecision(4) << 100.0 * xsize << " " << std::setw(7) << std::right << std::fixed
            << std::setprecision(4) << 100.0 * ysize << " "
            << "  0.2000 " << std::setw(6) << std::right << std::fixed << std::setprecision(2) << 100.0 * detd << " "
            << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.X() << " "
            << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.Y() << " "
            << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.Z() << " "
            << std::setw(8) << std::right << std::fixed << std::setprecision(5) << base.X() << " " << std::setw(8)
            << std::right << std::fixed << std::setprecision(5) << base.Y() << " " << std::setw(8) << std::right
            << std::fixed << std::setprecision(5) << base.Z() << " " << std::setw(8) << std::right << std::fixed
            << std::setprecision(5) << up.X() << " " << std::setw(8) << std::right << std::fixed << std::setprecision(5)
            << up.Y() << " " << std::setw(8) << std::right << std::fixed << std::setprecision(5) << up.Z() << " \n";
 
      } else
        g_log.warning() << "Information about detector module " << bankName << " not found and recognised\n";
    }
  }
  // HKL's are flipped by -1 because of the internal Q convention
  // unless Crystallography convention
  double qSign = -1.0;
  if (m_ws->getConvention() == "Crystallography")
    qSign = 1.0;
 
  // Save all Peaks
  // Go in order of run numbers
  int sequenceNumber = appendPeakNumb;
  for (const auto &runBankMap : runMap) {
    // Start of a new run
    const int run = runBankMap.first;
    const auto &bankMap = runBankMap.second;
 
    for (const auto &bankIDs : bankMap) {
      // Start of a new bank.
      const int bank = bankIDs.first;
      const auto &ids = bankIDs.second;
 
      if (!ids.empty()) {
        // Write the bank header
        out << "0  NRUN DETNUM     CHI      PHI    OMEGA       MONCNT\n";
        out << "1 " << std::setw(5) << run << std::setw(7) << std::right << bank;
 
        // Determine goniometer angles by calculating from the goniometer matrix
        // of a peak in the list
        Goniometer gon(peaks[ids[0]].getGoniometerMatrix());
        std::vector<double> angles = gon.getEulerAngles("yzy");
 
        double phi = angles[2];
        double chi = angles[1];
        double omega = angles[0];
 
        out << std::setw(8) << std::fixed << std::setprecision(2) << chi << " ";
        out << std::setw(8) << std::fixed << std::setprecision(2) << phi << " ";
        out << std::setw(8) << std::fixed << std::setprecision(2) << omega << " ";
 
        // Get the monitor count from the first peak (should all be the same for
        // one run)
        const size_t first_peak_index = ids[0];
        const auto &first_peak = peaks[first_peak_index];
        const double monct = first_peak.getMonitorCount();
        out << std::setw(12) << static_cast<int>(monct) << '\n';
        out << header << '\n';
 
        // Go through each peak at this run / bank
        for (auto wi : ids) {
          const auto &peak = peaks[wi];
 
          // Sequence (run) number
          std::string firstNumber = "3";
          if (m_isModulatedStructure) {
            firstNumber = "9";
          }
          if (renumber) {
            out << firstNumber << std::setw(7) << sequenceNumber;
            sequenceNumber++;
          } else {
            out << firstNumber << std::setw(7) << peak.getPeakNumber() + appendPeakNumb;
          }
 
          // HKL's are flipped by -1 because of the internal Q convention
          // unless Crystallography convention
          if (m_isModulatedStructure) {
            const V3D mod = peak.getIntMNP();
            const auto intHKL = peak.getIntHKL();
            out << std::setw(5) << Utils::round(qSign * intHKL.X()) << std::setw(5) << Utils::round(qSign * intHKL.Y())
                << std::setw(5) << Utils::round(qSign * intHKL.Z());
 
            out << std::setw(5) << Utils::round(qSign * mod[0]) << std::setw(5) << Utils::round(qSign * mod[1])
                << std::setw(5) << Utils::round(qSign * mod[2]);
          } else {
            out << std::setw(5) << Utils::round(qSign * peak.getH()) << std::setw(5)
                << Utils::round(qSign * peak.getK()) << std::setw(5) << Utils::round(qSign * peak.getL());
          }
 
          // Row/column
          out << std::setw(8) << std::fixed << std::setprecision(2) << static_cast<double>(peak.getCol()) << " ";
 
          out << std::setw(8) << std::fixed << std::setprecision(2) << static_cast<double>(peak.getRow()) << " ";
 
          out << std::setw(8) << std::fixed << std::setprecision(0) << peak.getTOF() << " ";
 
          out << std::setw(9) << std::fixed << std::setprecision(3) << (peak.getL2() * 100.0) << " ";
 
          // This is the scattered beam direction
          const V3D dir = peak.getDetPos() - inst->getSample()->getPos();
          double scattering, azimuth;
 
          // Two-theta = polar angle = scattering angle = between +Z vector and
          // the scattered beam
          scattering = dir.angle(V3D(0.0, 0.0, 1.0));
 
          // "Azimuthal" angle: project the scattered beam direction onto the XY
          // plane,
          // and calculate the angle between that and the +X axis (right-handed)
          azimuth = atan2(dir.Y(), dir.X());
 
          out << std::setw(9) << std::fixed << std::setprecision(5) << scattering << " "; // two-theta scattering
 
          out << std::setw(9) << std::fixed << std::setprecision(5) << azimuth << " ";
 
          out << std::setw(10) << std::fixed << std::setprecision(6) << peak.getWavelength() << " ";
 
          out << std::setw(9) << std::fixed << std::setprecision(4) << peak.getDSpacing() << " ";
 
          out << std::setw(8) << std::fixed << std::setprecision(0) << int(peak.getBinCount()) << " ";
 
          out << std::setw(10) << std::fixed << std::setprecision(2) << peak.getIntensity() << " ";
 
          out << std::setw(7) << std::fixed << std::setprecision(2) << peak.getSigmaIntensity() << " ";
 
          int thisReflag = 310;
          out << std::setw(5) << thisReflag;
 
          out << '\n';
 
          Workspace2D_sptr wsProfile2D = getProperty("ProfileWorkspace");
          if (wsProfile2D) {
            out << "8";
            const auto &yValues = wsProfile2D->y(wi);
            for (size_t j = 0; j < yValues.size(); j++) {
              out << std::setw(8) << static_cast<int>(yValues[j]);
              if ((j + 1) % 10 == 0) {
                out << '\n';
                if (j + 1 != yValues.size())
                  out << "8";
              }
            }
          }
        }
      }
    }
  }
 
  out.flush();
  out.close();
}
 
bool SaveIsawPeaks::bankMasked(size_t componentIndex, const Geometry::DetectorInfo &detectorInfo) {
  const auto &componentInfo = this->m_ws->componentInfo();
  auto children = componentInfo.children(componentIndex);
 
  if (!children.empty() && componentInfo.name(children[0]) == "sixteenpack") {
    children = componentInfo.children(children[0]);
  }
 
  for (const auto &colIndex : children) { // cppcheck-suppress useStlAlgorithm
    auto grandchildren = componentInfo.children(colIndex);
    for (const auto &rowIndex : grandchildren) { // cppcheck-suppress useStlAlgorithm
      if (componentInfo.isDetector(rowIndex)) {
        const auto detID = detectorInfo.detid(rowIndex);
        if (detID < 0)
          continue;
        const auto detIndex = detectorInfo.indexOf(detID);
        if (!detectorInfo.isMasked(detIndex)) {
          return false;
        }
      }
    }
  }
  return true;
}
 
V3D SaveIsawPeaks::findPixelPos(const std::string &bankName, int col, int row) {
  auto parent = inst->getComponentByName(bankName);
  if (parent->type() == "RectangularDetector") {
    const auto &componentInfo = this->m_ws->componentInfo();
    const size_t bankIndex = componentInfo.indexOf(parent->getComponentID());
    return componentInfo.position(componentInfo.detectorIndexAtXYZ(bankIndex, col, row, 0));
  } else {
    const auto &componentInfo = this->m_ws->componentInfo();
    const size_t parentIndex = componentInfo.indexOfAny(bankName);
    auto children = componentInfo.children(parentIndex);
 
    if (!children.empty() && componentInfo.name(children[0]) == "sixteenpack") {
      children = componentInfo.children(children[0]);
    }
    int col0 = col - 1;
    // WISH detectors are in bank in this order in instrument
    if (inst->getName() == "WISH")
      col0 = (col % 2 == 0 ? col / 2 + 75 : (col - 1) / 2);
 
    auto grandchildren = componentInfo.children(children[col0]);
    // NOTE: query ComponentInfo rather than dereferencing componentID(). The latter yields the *base* component,
    //       whose getPos() ignores the ParameterMap and so reports uncalibrated IDF geometry.
    return componentInfo.position(grandchildren[row - 1]);
  }
}
void SaveIsawPeaks::sizeBanks(const std::string &bankName, int &NCOLS, int &NROWS, double &xsize, double &ysize) {
  if (bankName == "None")
    return;
  const auto &componentInfo = this->m_ws->componentInfo();
  size_t parentIndex;
  try {
    parentIndex = componentInfo.indexOfAny(bankName);
  } catch (std::invalid_argument &) {
    return;
  }
  if (componentInfo.isGridDetector(parentIndex)) {
    const auto grid = componentInfo.pixelGridComponent(parentIndex);
 
    NCOLS = grid.nX;
    NROWS = grid.nY;
    xsize = grid.nX * grid.xStep;
    ysize = grid.nY * grid.yStep;
  } else {
    auto children = componentInfo.children(parentIndex);
 
    if (!children.empty() && componentInfo.name(children[0]) == "sixteenpack") {
      children = componentInfo.children(children[0]);
    }
 
    auto grandchildren = componentInfo.children(children[0]);
    NROWS = static_cast<int>(grandchildren.size());
    NCOLS = static_cast<int>(children.size());
    // NOTE: measure via ComponentInfo::position(), not by dereferencing componentID(). The latter yields the *base*
    //       component, whose getDistance() ignores the ParameterMap and so reports uncalibrated IDF geometry.
    xsize = componentInfo.position(children[0]).distance(componentInfo.position(children[NCOLS - 1]));
    ysize = componentInfo.position(grandchildren[0]).distance(componentInfo.position(grandchildren[NROWS - 1]));
  }
}
void SaveIsawPeaks::writeOffsets(std::ofstream &out, double qSign, const std::vector<double> &offset) {
  for (size_t i = 0; i < 3; i++) {
    out << std::setw(12) << std::fixed << std::setprecision(6) << qSign * offset[i] << " ";
  }
}
} // namespace Mantid::Crystal