Newer
Older
/*WIKI*
TODO: Enter a full wiki-markup description of your algorithm here. You can then use the Build/wiki_maker.py script to generate your full wiki page.
*WIKI*/
#include "MantidDataHandling/LoadPSI.h"
#include "MantidAPI/FileProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/LoadAlgorithmFactory.h"
#include "MantidAPI/Progress.h"
#include "MantidGeometry/Instrument.h"
#include <limits>
#include <algorithm>
#include <iostream>
#include <vector>
#include <cmath>
namespace Mantid {
namespace DataHandling {
using namespace Kernel;
using namespace API;
using namespace NeXus;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(LoadPSI)
//register the algorithm into loadalgorithm factory
DECLARE_LOADALGORITHM(LoadPSI)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
LoadPSI::LoadPSI() {
m_instrumentName = "";
supportedInstruments.push_back("FOCUS");
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
LoadPSI::~LoadPSI() {
}
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string LoadPSI::name() const {
return "LoadPSI";
}
;
/// Algorithm's version for identification. @see Algorithm::version
int LoadPSI::version() const {
return 1;
}
;
/// Algorithm's category for identification. @see Algorithm::category
const std::string LoadPSI::category() const {
return "DataHandling";
}
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void LoadPSI::initDocs() {
this->setWikiSummary("Loads PSI nexus file.");
this->setOptionalMessage("Loads PSI nexus file.");
}
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
bool LoadPSI::quickFileCheck(const std::string& filePath, size_t nread,
const file_header& header) {
std::string extn = extension(filePath);
bool bnexs(false);
(!extn.compare("nxs") || !extn.compare("nx5")) ? bnexs = true : bnexs =
false;
/*
* HDF files have magic cookie in the first 4 bytes
*/
if (((nread >= sizeof(unsigned))
&& (ntohl(header.four_bytes) == g_hdf_cookie)) || bnexs) {
//hdf
return true;
} else if ((nread >= sizeof(g_hdf5_signature))
&& (!memcmp(header.full_hdr, g_hdf5_signature,
sizeof(g_hdf5_signature)))) {
//hdf5
return true;
}
return false;
}
/**
* Checks the file by opening it and reading few lines
* @param filePath :: name of the file inluding its path
* @return an integer value how much this algorithm can load the file
*/
int LoadPSI::fileCheck(const std::string& filePath) {
// Create the root Nexus class
NXRoot root(filePath);
NXEntry entry = root.openFirstEntry();
std::string instrumentName = getInstrumentName(entry);
if (std::find(supportedInstruments.begin(), supportedInstruments.end(),
instrumentName) != supportedInstruments.end()) {
// FOUND
return 80;
}
return 0;
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void LoadPSI::init() {
std::vector<std::string> exts;
exts.push_back(".nxs");
exts.push_back(".hdf");
declareProperty(new FileProperty("Filename", "", FileProperty::Load, exts),
"The name of the Nexus file to load");
declareProperty(
new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
"The name to use for the output workspace");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void LoadPSI::exec() {
std::string filename = getPropertyValue("Filename");
NXRoot root(filename);
NXEntry entry = root.openFirstEntry();
loadExperimentDetails(entry);
runLoadInstrument();
setProperty("OutputWorkspace", m_localWorkspace);
}
void LoadPSI::setInstrumentName(NeXus::NXEntry& entry) {
m_instrumentName = getInstrumentName(entry);
if (m_instrumentName == "") {
std::string message(
"Cannot read the instrument name from the Nexus file!");
g_log.error(message);
throw std::runtime_error(message);
}
}
std::string LoadPSI::getInstrumentName(NeXus::NXEntry& entry) {
// format: /entry0/?????/name
std::string instrumentName = "";
std::vector<NXClassInfo> v = entry.groups();
for (auto it = v.begin(); it < v.end(); it++) {
if (it->nxclass == "NXinstrument") {
m_nexusInstrumentEntryName = it->nxname;
std::string insNamePath = m_nexusInstrumentEntryName + "/name";
instrumentName = entry.getString(insNamePath);
g_log.debug() << "Instrument Name: " << instrumentName
<< " in NxPath: " << insNamePath << std::endl;
break;
//std::replace( instrumentName.begin(), instrumentName.end(), ' ', '_'); // replace all ' ' to '_'
long unsigned int pos = instrumentName.find(" ");
instrumentName = instrumentName.substr (0,pos);
}
void LoadPSI::initWorkSpace(NeXus::NXEntry& entry) {
// read in the data
NXData dataGroup = entry.openNXData("merged");
NXInt data = dataGroup.openIntData();
m_numberOfTubes = static_cast<size_t>(data.dim0());
m_numberOfPixelsPerTube = 1;
m_numberOfChannels = static_cast<size_t>(data.dim1());
// dim0 * m_numberOfPixelsPerTube is the total number of detectors
m_numberOfHistograms = m_numberOfTubes * m_numberOfPixelsPerTube;
g_log.debug() << "NumberOfTubes: " << m_numberOfTubes << std::endl;
g_log.debug() << "NumberOfPixelsPerTube: " << m_numberOfPixelsPerTube
<< std::endl;
g_log.debug() << "NumberOfChannels: " << m_numberOfChannels << std::endl;
// Now create the output workspace
// Might need to get this value from the number of monitors in the Nexus file
// params:
// workspace type,
// total number of spectra + (number of monitors = 0),
// bin boundaries = m_numberOfChannels + 1
// Z/time dimension
m_localWorkspace = WorkspaceFactory::Instance().create("Workspace2D",
m_numberOfHistograms, m_numberOfChannels + 1, m_numberOfChannels);
m_localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create(
"TOF");
m_localWorkspace->setYUnitLabel("Counts");
}
void LoadPSI::loadDataIntoTheWorkSpace(NeXus::NXEntry& entry) {
// read in the data
NXData dataGroup = entry.openNXData("merged");
NXInt data = dataGroup.openIntData();
data.load();
NXFloat timeBinning = entry.openNXFloat("merged/time_binning");
timeBinning.load();
size_t numberOfBins = static_cast<size_t>(timeBinning.dim0()) + 1; // boundaries
g_log.debug() << "Number of bins: " << numberOfBins << std::endl;
float* timeBinning_p = &timeBinning[0];
std::vector<double> timeBinningTmp(numberOfBins);
timeBinningTmp.assign(timeBinning_p,timeBinning_p + numberOfBins);
timeBinningTmp[numberOfBins-1] = timeBinningTmp[numberOfBins-2]+timeBinningTmp[1]-timeBinningTmp[0];
m_localWorkspace->dataX(0).assign(timeBinningTmp.begin(),timeBinningTmp.end());
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
Progress progress(this, 0, 1, m_numberOfTubes * m_numberOfPixelsPerTube);
size_t spec = 0;
for (size_t i = 0; i < m_numberOfTubes; ++i) {
for (size_t j = 0; j < m_numberOfPixelsPerTube; ++j) {
if (spec > 0) {
// just copy the time binning axis to every spectra
m_localWorkspace->dataX(spec) = m_localWorkspace->readX(0);
}
// Assign Y
int* data_p = &data(static_cast<int>(i), static_cast<int>(j));
m_localWorkspace->dataY(spec).assign(data_p,
data_p + m_numberOfChannels);
// Assign Error
MantidVec& E = m_localWorkspace->dataE(spec);
std::transform(data_p, data_p + m_numberOfChannels, E.begin(),
LoadPSI::calculateError);
++spec;
progress.report();
}
}
g_log.debug() << "Data loading inti WS done...." << std::endl;
}
void LoadPSI::loadRunDetails(NXEntry & entry) {
API::Run & runDetails = m_localWorkspace->mutableRun();
// int runNum = entry.getInt("run_number");
// std::string run_num = boost::lexical_cast<std::string>(runNum);
// runDetails.addProperty("run_number", run_num);
std::string start_time = entry.getString("start_time");
//start_time = getDateTimeInIsoFormat(start_time);
runDetails.addProperty("run_start", start_time);
std::string end_time = entry.getString("end_time");
//end_time = getDateTimeInIsoFormat(end_time);
runDetails.addProperty("run_end", end_time);
double wavelength = entry.getFloat(m_nexusInstrumentEntryName + "/monochromator/lambda");
runDetails.addProperty<double>("wavelength", wavelength);
double energy = entry.getFloat(m_nexusInstrumentEntryName + "/monochromator/energy");
runDetails.addProperty<double>("Ei", energy,true); //overwrite
std::string title = entry.getString("title");
runDetails.addProperty("title", title);
m_localWorkspace->setTitle(title);
}
/*
* Load data about the Experiment.
*
* TODO: This is very incomplete. In ISIS they much more info in the nexus file than ILL.
*
* @param entry :: The Nexus entry
*/
void LoadPSI::loadExperimentDetails(NXEntry & entry) {
// TODO: Do the rest
// Pick out the geometry information
// std::string description = boost::lexical_cast<std::string>(
// entry.getFloat("sample/description"));
//
// m_localWorkspace->mutableSample().setName(description);
// m_localWorkspace->mutableSample().setThickness(static_cast<double> (isis_raw->spb.e_thick));
// m_localWorkspace->mutableSample().setHeight(static_cast<double> (isis_raw->spb.e_height));
// m_localWorkspace->mutableSample().setWidth(static_cast<double> (isis_raw->spb.e_width));
}
/**
* Run the Child Algorithm LoadInstrument.
*/
void LoadPSI::runLoadInstrument() {
IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
// Now execute the Child Algorithm. Catch and log any error, but don't stop.
try {
// TODO: depending on the m_numberOfPixelsPerTube we might need to load a different IDF
loadInst->setPropertyValue("InstrumentName", m_instrumentName);
loadInst->setProperty<MatrixWorkspace_sptr>("Workspace",
m_localWorkspace);
loadInst->execute();
} catch (...) {
g_log.information("Cannot load the instrument definition.");
}
}
} // namespace DataHandling
} // namespace Mantid