SaveIsawDetCal.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 "MantidDataHandling/SaveIsawDetCal.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Workspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/ArrayProperty.h"
 
#include "MantidKernel/Strings.h"
 
#include <boost/algorithm/string/trim.hpp>
#include <filesystem>
#include <fstream>
 
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using Mantid::Types::Core::DateAndTime;
using std::string;
 
namespace Mantid::DataHandling {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SaveIsawDetCal)
 
 
void SaveIsawDetCal::init() {
  declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input),
                  "An input workspace.");
 
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Save, ".DetCal"),
                  "Path to an ISAW-style .detcal file to save.");
 
  declareProperty("TimeOffset", 0.0, "Offsets to be applied to times");
  declareProperty(std::make_unique<ArrayProperty<string>>("BankNames", Direction::Input),
                  "Optional: Only select the specified banks");
  declareProperty("AppendFile", false,
                  "Append to file if true.\n"
                  "If false, new file (default).");
}
 
void SaveIsawDetCal::exec() {
  std::string filename = getPropertyValue("Filename");
  // MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
 
  Workspace_sptr ws1 = getProperty("InputWorkspace");
  ExperimentInfo_sptr ws = std::dynamic_pointer_cast<ExperimentInfo>(ws1);
 
  double T0 = getProperty("TimeOffset");
  const API::Run &run = ws->run();
  // Use T0 from workspace if T0 not specified in input
  if (T0 == 0.0 && run.hasProperty("T0")) {
    T0 = run.getPropertyValueAsType<double>("T0");
    if (T0 != 0) {
      g_log.notice() << "T0 = " << T0 << '\n';
    }
  }
 
  std::vector<std::string> bankNames = getProperty("BankNames");
 
  inst = ws->getInstrument();
  if (!inst)
    throw std::runtime_error("No instrument in the Workspace. Cannot save DetCal file.");
  // 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";
 
  std::set<int> uniqueBanks;
  if (bankNames.empty()) {
    // Get all components using ComponentInfo
    const auto &componentInfo = ws->componentInfo();
    for (size_t i = 0; i < componentInfo.size(); ++i) {
      std::string bankName = componentInfo.name(i);
      boost::trim(bankName);
      boost::erase_all(bankName, bankPart);
      int bank = 0;
      Strings::convert(bankName, bank);
      if (bank == 0)
        continue;
      // Track unique bank numbers
      uniqueBanks.insert(bank);
    }
  } else {
    for (auto bankName : bankNames) {
      boost::trim(bankName);
      boost::erase_all(bankName, bankPart);
      int bank = 0;
      Strings::convert(bankName, bank);
      if (bank == 0)
        continue;
      // Track unique bank numbers
      uniqueBanks.insert(bank);
    }
  }
 
  double l1;
  V3D beamline;
  double beamline_norm;
  V3D samplePos;
  inst->getInstrumentParameters(l1, beamline, beamline_norm, samplePos);
 
  std::ofstream out;
  bool append = getProperty("AppendFile");
 
  // do not append if file does not exist
  if (!std::filesystem::exists(std::filesystem::path(filename.c_str())))
    append = false;
 
  if (append) {
    out.open(filename.c_str(), std::ios::app);
  } else {
    out.open(filename.c_str());
    out << "# NEW CALIBRATION FILE FORMAT (in NeXus/SNS coordinates):\n";
    out << "# Lengths are in centimeters.\n";
    out << "# Base and up give directions of unit vectors for a local \n";
    out << "# x,y coordinate system on the face of the detector.\n";
    out << "#\n";
    out << "#\n";
    out << "# " << DateAndTime::getCurrentTime().toISO8601String() << '\n';
 
    out << "6         L1     T0_SHIFT\n";
    out << "7 " << std::setw(10);
    out << std::setprecision(4) << std::fixed << (l1 * 100);
    out << std::setw(13) << std::setprecision(3) << T0 << '\n';
 
    out << "4 DETNUM  NROWS  NCOLS   WIDTH   HEIGHT   DEPTH   DETD   CenterX "
           "  CenterY   CenterZ    BaseX    BaseY    BaseZ      UpX      UpY "
           "     UpZ\n";
  }
  // Here would save each detector...
  const auto &componentInfo = ws->componentInfo();
  std::set<int>::iterator it;
  for (it = uniqueBanks.begin(); it != uniqueBanks.end(); ++it) {
    // Build up the bank name
    int bank = *it;
    std::ostringstream mess;
    if (bankPart == "WISHpanel" && bank < 10)
      mess << bankPart << "0" << bank;
    else
      mess << bankPart << 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 = IComponent_const_sptr(componentInfo.componentID(children[0]), NoDeleting());
      }
    }
    if (det) {
      // Center of the detector
      V3D center = det->getPos();
 
      // Distance to center of detector
      double detd = (center - inst->getSample()->getPos()).norm();
      int NCOLS, NROWS;
      double xsize, ysize;
      sizeBanks(bankName, NCOLS, NROWS, xsize, ysize, componentInfo);
      // Base unit vector (along the horizontal, X axis)
      // Since WISH banks are curved use center and increment 2 for tubedown
      int midX = NCOLS / 2;
      int midY = NROWS / 2;
      V3D base =
          findPixelPos(bankName, midX + 2, midY, componentInfo) - findPixelPos(bankName, midX, midY, componentInfo);
      base.normalize();
 
      // Up unit vector (along the vertical, Y axis)
      V3D up =
          findPixelPos(bankName, midX, midY + 1, componentInfo) - findPixelPos(bankName, midX, midY, componentInfo);
      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";
  }
 
  out.close();
}
 
V3D SaveIsawDetCal::findPixelPos(const std::string &bankName, int col, int row, const ComponentInfo &componentInfo) {
  std::shared_ptr<const IComponent> parent = inst->getComponentByName(bankName);
  if (parent->type() == "RectangularDetector") {
    std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
 
    std::shared_ptr<Detector> pixel = RDet->getAtXY(col, row);
    return pixel->getPos();
  } else {
    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);
    const auto grandchildren = componentInfo.children(children[col0]);
    return componentInfo.position(grandchildren[row - 1]);
  }
}
 
void SaveIsawDetCal::sizeBanks(const std::string &bankName, int &NCOLS, int &NROWS, double &xsize, double &ysize,
                               const ComponentInfo &componentInfo) {
  if (bankName == "None")
    return;
  std::shared_ptr<const IComponent> parent = inst->getComponentByName(bankName);
  if (parent->type() == "RectangularDetector") {
    std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
 
    NCOLS = RDet->xpixels();
    NROWS = RDet->ypixels();
    xsize = RDet->xsize();
    ysize = RDet->ysize();
  } else {
    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]);
    }
    const auto grandchildren = componentInfo.children(children[0]);
    NROWS = static_cast<int>(grandchildren.size());
    NCOLS = static_cast<int>(children.size());
    const auto &first = componentInfo.componentID(children[0]);
    const auto &last = componentInfo.componentID(children[NCOLS - 1]);
    xsize = first->getDistance(*last);
    const auto &firstGrand = componentInfo.componentID(grandchildren[0]);
    const auto &lastGrand = componentInfo.componentID(grandchildren[NROWS - 1]);
    ysize = firstGrand->getDistance(*lastGrand);
  }
}
} // namespace Mantid::DataHandling