Skip to content
Snippets Groups Projects
Commit 2eb7b1c6 authored by Lynch, Vickie's avatar Lynch, Vickie
Browse files

Merge pull request #474 from mantidproject/feature/11288_refactor_convertspice

Refactor ConvertSpiceDataToRealWorkspace
parents e1a72e59 c4141a12
No related branches found
No related tags found
No related merge requests found
......@@ -10,7 +10,8 @@
namespace Mantid {
namespace MDAlgorithms {
/** LoadHFIRPDD : TODO: DESCRIPTION
/** ConvertSpiceDataToRealSpace : Convert data from SPICE file to singals
in real space contained in MDEventWrokspaces
Copyright © 2014 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
......@@ -60,6 +61,9 @@ public:
// virtual int confidence(Kernel::FileDescriptor &descriptor) const;
private:
/// Typdef for the white-space separated file data type.
typedef std::deque<std::string> DataCollectionType;
/// Initialisation code
void init();
......@@ -70,17 +74,21 @@ private:
DataObjects::TableWorkspace_sptr
loadSpiceData(const std::string &spicefilename);
/// Convert to MD workspaces
/// Parse data table workspace to a vector of matrix workspaces
std::vector<API::MatrixWorkspace_sptr> convertToMatrixWorkspace(
DataObjects::TableWorkspace_sptr tablews,
API::MatrixWorkspace_const_sptr parentws, Kernel::DateAndTime runstart,
std::map<std::string, std::vector<double> > &logvecmap,
std::vector<Kernel::DateAndTime> &vectimes);
/// Create an MDEventWorspace by converting vector of matrix workspace data
API::IMDEventWorkspace_sptr
convertToMDEventWS(const std::vector<API::MatrixWorkspace_sptr> &vec_ws2d);
createDataMDWorkspace(const std::vector<API::MatrixWorkspace_sptr> &vec_ws2d);
/// Parse data table workspace to a vector of matrix workspaces
std::vector<API::MatrixWorkspace_sptr>
convertToWorkspaces(DataObjects::TableWorkspace_sptr tablews,
API::MatrixWorkspace_const_sptr parentws,
Kernel::DateAndTime runstart,
std::map<std::string, std::vector<double> > &logvecmap,
std::vector<Kernel::DateAndTime> &vectimes);
/// Create an MDWorkspace for monitor counts
API::IMDEventWorkspace_sptr createMonitorMDWorkspace(
const std::vector<API::MatrixWorkspace_sptr> vec_ws2d,
const std::vector<double> &vecmonitor);
/// Read parameter information from table workspace
void readTableInfo(DataObjects::TableWorkspace_const_sptr tablews,
......@@ -113,16 +121,20 @@ private:
const std::map<std::string, std::vector<double> > &logvecmap,
const std::vector<Kernel::DateAndTime> &vectimes);
/// Create an MDWorkspace for monitor counts
API::IMDEventWorkspace_sptr createMonitorMDWorkspace(
const std::vector<API::MatrixWorkspace_sptr> vec_ws2d,
const std::vector<double> &vecmonitor);
/// Name of instrument
std::string m_instrumentName;
/// Number of detectors
size_t m_numSpec;
/// x-y-z-value minimum
std::vector<double> m_extentMins;
/// x-y-z value maximum
std::vector<double> m_extentMaxs;
/// Number of bins
std::vector<size_t> m_numBins;
/// Dimension of the output MDEventWorkspace
size_t m_nDimensions;
};
} // namespace DataHandling
......
......@@ -7,6 +7,12 @@
#include "MantidGeometry/IDetector.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDEventInserter.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidDataObjects/MDEvent.h"
#include <boost/algorithm/string/predicate.hpp>
#include <Poco/TemporaryFile.h>
......@@ -17,6 +23,8 @@ namespace MDAlgorithms {
using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
DECLARE_ALGORITHM(ConvertSpiceDataToRealSpace)
......@@ -24,14 +32,14 @@ DECLARE_ALGORITHM(ConvertSpiceDataToRealSpace)
/** Constructor
*/
ConvertSpiceDataToRealSpace::ConvertSpiceDataToRealSpace()
: m_instrumentName(""), m_numSpec(0) {}
: m_instrumentName(""), m_numSpec(0), m_nDimensions(3) {}
//------------------------------------------------------------------------------------------------
/** Destructor
*/
ConvertSpiceDataToRealSpace::~ConvertSpiceDataToRealSpace() {}
//----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Init
*/
void ConvertSpiceDataToRealSpace::init() {
......@@ -84,11 +92,10 @@ void ConvertSpiceDataToRealSpace::init() {
"Name to use for the output workspace.");
}
//----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Exec
*/
void ConvertSpiceDataToRealSpace::exec() {
// Process inputs
DataObjects::TableWorkspace_sptr dataTableWS = getProperty("InputWorkspace");
MatrixWorkspace_const_sptr parentWS = getProperty("RunInfoWorkspace");
......@@ -138,12 +145,35 @@ void ConvertSpiceDataToRealSpace::exec() {
// Convert table workspace to a list of 2D workspaces
std::map<std::string, std::vector<double> > logvecmap;
std::vector<Kernel::DateAndTime> vectimes;
std::vector<MatrixWorkspace_sptr> vec_ws2d =
convertToWorkspaces(dataTableWS, parentWS, runstart, logvecmap, vectimes);
// Set up range for x/y/z
m_extentMins.resize(3);
m_extentMaxs.resize(3);
for (size_t i = 0; i < 3; ++i) {
m_extentMins[i] = DBL_MAX;
m_extentMaxs[i] = -DBL_MAX;
}
std::vector<MatrixWorkspace_sptr> vec_ws2d = convertToMatrixWorkspace(
dataTableWS, parentWS, runstart, logvecmap, vectimes);
// check range for x/y/z
m_numBins.resize(3);
for (size_t d = 0; d < 3; ++d) {
if (fabs(m_extentMins[d] - m_extentMaxs[d]) < 1.0E-6) {
// Range is too small so treat it as 1 value
double mvalue = m_extentMins[d];
m_extentMins[d] = mvalue - 0.1;
m_extentMaxs[d] = mvalue + 0.1;
m_numBins[d] = 1;
} else {
m_numBins[d] = 100;
}
}
// Convert to MD workspaces
g_log.debug("About to converting to workspaces done!");
IMDEventWorkspace_sptr m_mdEventWS = convertToMDEventWS(vec_ws2d);
IMDEventWorkspace_sptr m_mdEventWS = createDataMDWorkspace(vec_ws2d);
std::string monitorlogname = getProperty("MonitorCountsLogName");
IMDEventWorkspace_sptr mdMonitorWS =
createMonitorMDWorkspace(vec_ws2d, logvecmap[monitorlogname]);
......@@ -159,7 +189,7 @@ void ConvertSpiceDataToRealSpace::exec() {
setProperty("OutputMonitorWorkspace", mdMonitorWS);
}
//----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Convert runs/pts from table workspace to a list of workspace 2D
* @brief ConvertSpiceDataToRealSpace::convertToWorkspaces
* @param tablews
......@@ -170,7 +200,7 @@ void ConvertSpiceDataToRealSpace::exec() {
* @return
*/
std::vector<MatrixWorkspace_sptr>
ConvertSpiceDataToRealSpace::convertToWorkspaces(
ConvertSpiceDataToRealSpace::convertToMatrixWorkspace(
DataObjects::TableWorkspace_sptr tablews,
API::MatrixWorkspace_const_sptr parentws, Kernel::DateAndTime runstart,
std::map<std::string, std::vector<double> > &logvecmap,
......@@ -222,7 +252,7 @@ void ConvertSpiceDataToRealSpace::parseSampleLogs(
std::string logname = indexiter->first;
size_t icol = indexiter->second;
g_log.information() << " Parsing log " << logname << "\n";
g_log.debug() << " Parsing log " << logname << "\n";
std::vector<double> logvec(numrows);
for (size_t ir = 0; ir < numrows; ++ir) {
......@@ -236,7 +266,7 @@ void ConvertSpiceDataToRealSpace::parseSampleLogs(
return;
}
//----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Load one run of data to a new workspace
* @brief ConvertSpiceDataToRealSpace::loadRunToMatrixWS
* @param tablews
......@@ -296,16 +326,28 @@ MatrixWorkspace_sptr ConvertSpiceDataToRealSpace::loadRunToMatrixWS(
tempws = instloader->getProperty("Workspace");
// Import data
std::vector<double> pos(3);
for (size_t i = 0; i < m_numSpec; ++i) {
// get detector
Geometry::IDetector_const_sptr tmpdet = tempws->getDetector(i);
tempws->dataX(i)[0] = tmpdet->getPos().X();
tempws->dataX(i)[0] = tmpdet->getPos().X() + 0.01;
pos[0] = tmpdet->getPos().X();
pos[1] = tmpdet->getPos().Y();
pos[2] = tmpdet->getPos().Z();
tempws->dataX(i)[0] = pos[0];
tempws->dataX(i)[0] = pos[0] + 0.01;
double yvalue = tablews->cell<double>(irow, anodelist[i].second);
tempws->dataY(i)[0] = yvalue;
if (yvalue >= 1)
tempws->dataE(i)[0] = sqrt(yvalue);
else
tempws->dataE(i)[0] = 1;
// update X-range, Y-range and Z-range
for (size_t d = 0; d < 3; ++d) {
if (pos[d] < m_extentMins[d])
m_extentMins[d] = pos[d];
if (pos[d] > m_extentMaxs[d])
m_extentMaxs[d] = pos[d];
}
}
// Return duration
......@@ -314,7 +356,7 @@ MatrixWorkspace_sptr ConvertSpiceDataToRealSpace::loadRunToMatrixWS(
return tempws;
}
//----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Read table workspace's column information
* @brief ConvertSpiceDataToRealSpace::readTableInfo
* @param tablews
......@@ -394,175 +436,7 @@ void ConvertSpiceDataToRealSpace::readTableInfo(
return;
}
//----------------------------------------------------------------------------------------------
/** Convert to MD Event workspace
* @brief ConvertSpiceDataToRealSpace::convertToMDEventWS
* @param vec_ws2d
* @return
*/
IMDEventWorkspace_sptr ConvertSpiceDataToRealSpace::convertToMDEventWS(
const std::vector<MatrixWorkspace_sptr> &vec_ws2d) {
// Write the lsit of workspacs to a file to be loaded to an MD workspace
Poco::TemporaryFile tmpFile;
std::string tempFileName = tmpFile.path();
g_log.debug() << "Creating temporary MD Event file = " << tempFileName
<< "\n";
// Construct a file
std::ofstream myfile;
myfile.open(tempFileName.c_str());
myfile << "DIMENSIONS" << std::endl;
myfile << "x X m 100" << std::endl;
myfile << "y Y m 100" << std::endl;
myfile << "z Z m 100" << std::endl;
myfile << "# Signal, Error, RunId, DetectorId, coord1, coord2, ... to end of "
"coords" << std::endl;
myfile << "MDEVENTS" << std::endl;
if (vec_ws2d.size() > 0) {
Progress progress(this, 0, 1, vec_ws2d.size());
size_t detindex = 0;
for (auto it = vec_ws2d.begin(); it < vec_ws2d.end(); ++it) {
API::MatrixWorkspace_sptr thisWorkspace = *it;
int runnumber =
atoi(thisWorkspace->run().getProperty("run_number")->value().c_str());
std::size_t nHist = thisWorkspace->getNumberHistograms();
for (std::size_t i = 0; i < nHist; ++i) {
Geometry::IDetector_const_sptr det = thisWorkspace->getDetector(i);
const MantidVec &signal = thisWorkspace->readY(i);
const MantidVec &error = thisWorkspace->readE(i);
myfile << signal[0] << " ";
myfile << error[0] << " ";
myfile << runnumber << " ";
myfile << det->getID() + detindex << " ";
Kernel::V3D detPos = det->getPos();
myfile << detPos.X() << " ";
myfile << detPos.Y() << " ";
myfile << detPos.Z() << " ";
myfile << std::endl;
}
// Increment on detector IDs
detindex += nHist;
progress.report("Creating MD WS");
}
myfile.close();
} else {
throw std::runtime_error(
"There is no MatrixWorkspace to construct MDWorkspace.");
}
// Import to MD Workspace
IAlgorithm_sptr importMDEWS = createChildAlgorithm("ImportMDEventWorkspace");
// Now execute the Child Algorithm.
try {
importMDEWS->setPropertyValue("Filename", tempFileName);
importMDEWS->setProperty("OutputWorkspace", "Test");
importMDEWS->executeAsChildAlg();
}
catch (std::exception &exc) {
throw std::runtime_error(
std::string("Error running ImportMDEventWorkspace: ") + exc.what());
}
IMDEventWorkspace_sptr workspace =
importMDEWS->getProperty("OutputWorkspace");
if (!workspace)
throw(std::runtime_error("Can not retrieve results of child algorithm "
"ImportMDEventWorkspace"));
return workspace;
}
//-----------------------------------------------------------------------------------------------
/** Create an MDWorkspace for monitoring counts.
* @brief LoadHFIRPDD::createMonitorMDWorkspace
* @param vec_ws2d
* @param vecmonitor
* @return
*/
IMDEventWorkspace_sptr ConvertSpiceDataToRealSpace::createMonitorMDWorkspace(
const std::vector<MatrixWorkspace_sptr> vec_ws2d,
const std::vector<double> &vecmonitor) {
// Write the lsit of workspacs to a file to be loaded to an MD workspace
Poco::TemporaryFile tmpFile;
std::string tempFileName = tmpFile.path();
g_log.debug() << "Creating temporary MD Event file for monitor counts = "
<< tempFileName << "\n";
// Construct a file
std::ofstream myfile;
myfile.open(tempFileName.c_str());
myfile << "DIMENSIONS" << std::endl;
myfile << "x X m 100" << std::endl;
myfile << "y Y m 100" << std::endl;
myfile << "z Z m 100" << std::endl;
myfile << "# Signal, Error, RunId, coord1, DetectorId, coord2, ... to end of "
"coords" << std::endl;
myfile << "MDEVENTS" << std::endl;
if (vec_ws2d.size() > 0) {
Progress progress(this, 0, 1, vec_ws2d.size());
size_t detindex = 0;
for (auto it = vec_ws2d.begin(); it < vec_ws2d.end(); ++it) {
API::MatrixWorkspace_sptr thisWorkspace = *it;
int runnumber =
atoi(thisWorkspace->run().getProperty("run_number")->value().c_str());
double signal = vecmonitor[static_cast<size_t>(it - vec_ws2d.begin())];
std::size_t nHist = thisWorkspace->getNumberHistograms();
for (std::size_t i = 0; i < nHist; ++i) {
Geometry::IDetector_const_sptr det = thisWorkspace->getDetector(i);
// const MantidVec &signal = thisWorkspace->readY(i);
const MantidVec &error = thisWorkspace->readE(i);
myfile << signal << " ";
myfile << error[0] << " ";
myfile << runnumber << " ";
myfile << det->getID() + detindex << " ";
Kernel::V3D detPos = det->getPos();
myfile << detPos.X() << " ";
myfile << detPos.Y() << " ";
myfile << detPos.Z() << " ";
myfile << std::endl;
}
// Increment on detector IDs
detindex += nHist;
progress.report("Creating MD WS");
}
myfile.close();
} else {
throw std::runtime_error(
"There is no MatrixWorkspace to construct MDWorkspace.");
}
// Import to MD Workspace
IAlgorithm_sptr importMDEWS = createChildAlgorithm("ImportMDEventWorkspace");
// Now execute the Child Algorithm.
try {
importMDEWS->setPropertyValue("Filename", tempFileName);
importMDEWS->setProperty("OutputWorkspace", "Test");
importMDEWS->executeAsChildAlg();
}
catch (std::exception &exc) {
throw std::runtime_error(
std::string("Error running ImportMDEventWorkspace: ") + exc.what());
}
IMDEventWorkspace_sptr workspace =
importMDEWS->getProperty("OutputWorkspace");
if (!workspace)
throw(std::runtime_error("Can not retrieve results of child algorithm "
"ImportMDEventWorkspace"));
return workspace;
}
//-----------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Create sample logs for MD workspace
* @brief LoadHFIRPDD::appendSampleLogs
* @param mdws
......@@ -650,7 +524,7 @@ void ConvertSpiceDataToRealSpace::appendSampleLogs(
return;
}
//---------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Add Experiment Info to the MDWorkspace. Add 1+N ExperimentInfo
* @brief ConvertSpiceDataToRealSpace::addExperimentInfos
* @param mdws
......@@ -683,5 +557,164 @@ void ConvertSpiceDataToRealSpace::addExperimentInfos(
return;
}
//------------------------------------------------------------------------------------------------
/** Convert to MD Event workspace
* @brief ConvertSpiceDataToRealSpace::convertToMDEventWS
* @param vec_ws2d
* @return
*/
IMDEventWorkspace_sptr ConvertSpiceDataToRealSpace::createDataMDWorkspace(
const std::vector<MatrixWorkspace_sptr> &vec_ws2d) {
// Create a target output workspace.
IMDEventWorkspace_sptr outWs =
MDEventFactory::CreateMDWorkspace(m_nDimensions, "MDEvent");
// Extract Dimensions and add to the output workspace.
std::vector<std::string> vec_ID(3);
vec_ID[0] = "x";
vec_ID[1] = "y";
vec_ID[2] = "z";
std::vector<std::string> vec_name(3);
vec_name[0] = "X";
vec_name[1] = "Y";
vec_name[2] = "Z";
// Add dimensions
for (size_t i = 0; i < m_nDimensions; ++i) {
std::string id = vec_ID[i];
std::string name = vec_name[i];
std::string units = "m";
// int nbins = 100;
for (size_t d = 0; d < 3; ++d)
g_log.debug() << "Direction " << d << ", Range = " << m_extentMins[d]
<< ", " << m_extentMaxs[d] << "\n";
outWs->addDimension(
Geometry::MDHistoDimension_sptr(new Geometry::MDHistoDimension(
id, name, units, static_cast<coord_t>(m_extentMins[i]),
static_cast<coord_t>(m_extentMaxs[i]), m_numBins[i])));
}
// Add events
// Creates a new instance of the MDEventInserter.
MDEventWorkspace<MDEvent<3>, 3>::sptr MDEW_MDEVENT_3 =
boost::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3> >(outWs);
MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(
MDEW_MDEVENT_3);
for (size_t iws = 0; iws < vec_ws2d.size(); ++iws) {
API::MatrixWorkspace_sptr thisWorkspace = vec_ws2d[iws];
short unsigned int runnumber = static_cast<short unsigned int>(
atoi(thisWorkspace->run().getProperty("run_number")->value().c_str()));
detid_t detindex = 0;
size_t nHist = thisWorkspace->getNumberHistograms();
for (std::size_t i = 0; i < nHist; ++i) {
// For each spectrum/detector
Geometry::IDetector_const_sptr det = thisWorkspace->getDetector(i);
const MantidVec &vecsignal = thisWorkspace->readY(i);
const MantidVec &vecerror = thisWorkspace->readE(i);
float signal = static_cast<float>(vecsignal[0]);
float error = static_cast<float>(vecerror[0]);
detid_t detid = det->getID() + detindex;
Kernel::V3D detPos = det->getPos();
double x = detPos.X();
double y = detPos.Y();
double z = detPos.Z();
std::vector<Mantid::coord_t> data(3);
data[0] = static_cast<float>(x);
data[1] = static_cast<float>(y);
data[2] = static_cast<float>(z);
inserter.insertMDEvent(signal, error * error, runnumber, detid,
data.data());
} // ENDFOR(spectrum)
} // ENDFOR (workspace)
return outWs;
}
//------------------------------------------------------------------------------------------------
/** Create an MDWorkspace for monitoring counts.
* @brief LoadHFIRPDD::createMonitorMDWorkspace
* @param vec_ws2d
* @param vecmonitor
* @return
*/
IMDEventWorkspace_sptr ConvertSpiceDataToRealSpace::createMonitorMDWorkspace(
const std::vector<MatrixWorkspace_sptr> vec_ws2d,
const std::vector<double> &vecmonitor) {
// Create a target output workspace.
IMDEventWorkspace_sptr outWs =
MDEventFactory::CreateMDWorkspace(m_nDimensions, "MDEvent");
// Extract Dimensions and add to the output workspace.
std::vector<std::string> vec_ID(3);
vec_ID[0] = "x";
vec_ID[1] = "y";
vec_ID[2] = "z";
std::vector<std::string> vec_name(3);
vec_name[0] = "X";
vec_name[1] = "Y";
vec_name[2] = "Z";
// Add dimensions
for (size_t i = 0; i < m_nDimensions; ++i) {
std::string id = vec_ID[i];
std::string name = vec_name[i];
std::string units = "m";
outWs->addDimension(
Geometry::MDHistoDimension_sptr(new Geometry::MDHistoDimension(
id, name, units, static_cast<coord_t>(m_extentMins[i]),
static_cast<coord_t>(m_extentMaxs[i]), m_numBins[i])));
}
// Add events
// Creates a new instance of the MDEventInserter.
MDEventWorkspace<MDEvent<3>, 3>::sptr MDEW_MDEVENT_3 =
boost::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3> >(outWs);
MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(
MDEW_MDEVENT_3);
for (size_t iws = 0; iws < vec_ws2d.size(); ++iws) {
API::MatrixWorkspace_sptr thisWorkspace = vec_ws2d[iws];
short unsigned int runnumber = static_cast<short unsigned int>(
atoi(thisWorkspace->run().getProperty("run_number")->value().c_str()));
detid_t detindex = 0;
float signal = static_cast<float>(vecmonitor[iws]);
float error = 1;
if (signal > 1)
error = std::sqrt(signal);
size_t nHist = thisWorkspace->getNumberHistograms();
for (std::size_t i = 0; i < nHist; ++i) {
// For each spectrum/detector
Geometry::IDetector_const_sptr det = thisWorkspace->getDetector(i);
detid_t detid = det->getID() + detindex;
Kernel::V3D detPos = det->getPos();
double x = detPos.X();
double y = detPos.Y();
double z = detPos.Z();
std::vector<Mantid::coord_t> data(3);
data[0] = static_cast<float>(x);
data[1] = static_cast<float>(y);
data[2] = static_cast<float>(z);
inserter.insertMDEvent(signal, error * error, runnumber, detid,
data.data());
} // ENDFOR(spectrum)
} // ENDFOR (workspace)
return outWs;
}
} // namespace DataHandling
} // namespace Mantid
......@@ -149,9 +149,9 @@ public:
Mantid::detid_t detid43 = mditer->getInnerDetectorID(43);
TS_ASSERT_EQUALS(detid43, 44);
Mantid::detid_t detid44 = mditer->getInnerDetectorID(44);
TS_ASSERT_EQUALS(detid44, 45);
TS_ASSERT_EQUALS(detid44, 1);
Mantid::detid_t detid61 = mditer->getInnerDetectorID(61);
TS_ASSERT_EQUALS(detid61, 62);
TS_ASSERT_EQUALS(detid61, 18);
// Run index
uint16_t run0 = mditer->getInnerRunIndex(0);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment