Skip to content

Commit

Permalink
Refs #5146 Refactor code into a few new methods
Browse files Browse the repository at this point in the history 
modified:   ../Framework/DataHandling/inc/MantidDataHandling/LoadSassena.h
modified:   ../Framework/DataHandling/src/LoadSassena.cpp
  • Loading branch information
@jmborr
jmborr committed May 22, 2012
1 parent 7ff0ca0 commit 639561c
Show file tree
Hide file tree
Showing 2 changed files with 195 additions and 101 deletions.
11 changes: 11 additions & 0 deletions Code/Mantid/Framework/DataHandling/inc/MantidDataHandling/LoadSassena.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include "MantidAPI/IDataFileChecker.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/Workspace2D.h"
#include <hdf5.h>

namespace Mantid
{
Expand Down Expand Up @@ -86,6 +87,16 @@ namespace Mantid
protected:
/// Add a workspace to the group and register in the analysis data service
void registerWorkspace(API::WorkspaceGroup_sptr gws, const std::string wsName, DataObjects::Workspace2D_sptr ws, const std::string &description);
/// Read info about one HDF5 dataset, log if error
void dataSetInfo( const hid_t& h5file, const std::string setName, hsize_t* dims);
/// Read dataset data to a buffer ot type double
void dataSetDouble( const hid_t& h5file, const std::string setName, double *buf );
/// Load qvectors dataset, calculate modulus of vectors
const MantidVec loadQvectors(const hid_t& h5file, API::WorkspaceGroup_sptr gws);
/// Load structure factor asa function of q-vector modulus
void loadFQ(const hid_t& h5file, API::WorkspaceGroup_sptr gws, const std::string setName, const MantidVec &qvmod);
/// Load time-dependent structure factor
void loadFQT(const hid_t& h5file, API::WorkspaceGroup_sptr gws, const std::string setName, const MantidVec &qvmod);

private:
/// Sets documentation strings for this algorithm
Expand Down
285 changes: 184 additions & 101 deletions Code/Mantid/Framework/DataHandling/src/LoadSassena.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@ It will create a workspace for each scattering intensity and one workspace for t
#include "MantidAPI/LoadAlgorithmFactory.h"
#include "MantidAPI/FileProperty.h"
#include "MantidDataObjects/Workspace2D.h"
#include <hdf5.h>
#include "MantidKernel/Exception.h"
//#include <hdf5.h>
#include <hdf5_hl.h>

namespace Mantid
Expand Down Expand Up @@ -58,16 +59,178 @@ int LoadSassena::fileCheck(const std::string &filePath)
{
int confidence(0);
hid_t h5file = H5Fopen(filePath.c_str(),H5F_ACC_RDONLY,H5P_DEFAULT);
if (H5LTfind_attribute(h5file,"sassena_version")==1) confidence = 99;
if (H5LTfind_attribute(h5file,"sassena_version")==1)
confidence = 99;
//else
//log.debug("no version attribute found");
return confidence;
}

/**
* Register a workspace in the Analysis Data Service and add it to the groupWorkspace
* @param gws pointer to WorkspaceGroup being filled
* @param wsName name of workspace to be added and registered
* @param ws pointer to workspace to be added and registered
* @param description string summarizing workspace contents
*/
void LoadSassena::registerWorkspace( API::WorkspaceGroup_sptr gws, const std::string wsName, DataObjects::Workspace2D_sptr ws, const std::string &description )
{
this->declareProperty(new API::WorkspaceProperty<DataObjects::Workspace2D>(wsName,"",Kernel::Direction::Output), description);
this->declareProperty(new API::WorkspaceProperty<DataObjects::Workspace2D>(wsName,wsName,Kernel::Direction::Output), description);
this->setProperty(wsName,ws);
gws->add(wsName); // Add the named workspace to the group
}

/**
* Read dataset dimensionality
* @param h5file file identifier
* @param setName string name of dataset
* @param buffer storing dimensionality
*/
void LoadSassena::dataSetInfo( const hid_t& h5file, const std::string setName, hsize_t* dims)
{
H5T_class_t class_id;
size_t type_size;
if( H5LTget_dataset_info( h5file, setName.c_str(), dims, &class_id, &type_size ) < 0 )
{
this->g_log.error("Unable to read "+setName+" dataset info");
throw Kernel::Exception::FileError("Unable to read "+setName+" dataset info:" , m_filename);
}
}

/**
* Read the dataset
* @param h5file file identifier
* @param setName string name of dataset
* @param buffer storing dataset
*/
void LoadSassena::dataSetDouble( const hid_t& h5file, const std::string setName, double *buf )
{
if( H5LTread_dataset_double(h5file,setName.c_str(),buf) < 0 )
{
this->g_log.error("Cannot read "+setName+" dataset");
throw Kernel::Exception::FileError("Unable to read "+setName+" dataset:" , m_filename);
}
}

/**
* load vectors onto a Workspace2D with 3 bins (the three components of the vectors)
* dataX for the origin of the vector (assumed (0,0,0) )
* dataY for the tip of the vector
* dataE is assumed (0,0,0), no errors
* @param h5file file identifier
* @param gws pointer to WorkspaceGroup being filled
*/
const MantidVec LoadSassena::loadQvectors(const hid_t& h5file, API::WorkspaceGroup_sptr gws)
{
const std::string gwsName = this->getPropertyValue("OutputWorkspace");
std::string setName("qvectors");

hsize_t dims[3];
this->dataSetInfo(h5file, setName, dims);
int nq = static_cast<int>( dims[0] ); //number of q-vectors
double* buf = new double[nq*3];
this->dataSetDouble(h5file,"qvectors",buf);

DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, 3, 3));
std::string wsName = gwsName + std::string("_") + setName;
ws->setTitle(wsName);

MantidVec qvmod; //store the modulus of the vector
double* curr = buf;
for(int iq=0; iq<nq; iq++){
MantidVec& Y = ws->dataY(iq);
Y.assign(curr,curr+3);
qvmod.push_back( sqrt( curr[0]*curr[0] + curr[1]*curr[1] + curr[2]*curr[2] ) );
curr += 3;
}
delete buf;
this->registerWorkspace(gws,wsName,ws, "X-axis: origin of Q-vectors; Y-axis: tip of Q-vectors");
return qvmod;
}

/**
* Create workspace to store the structure factor.
* First spectrum is the real part, second spectrum is the imaginary part
* X values are the modulus of the Q-vectors
* @param h5file file identifier
* @param gws pointer to WorkspaceGroup being filled
* @param setName string name of dataset
* @param qvmod vector of Q-vectors' moduli
*/
void LoadSassena::loadFQ(const hid_t& h5file, API::WorkspaceGroup_sptr gws, const std::string setName, const MantidVec &qvmod){
const std::string gwsName = this->getPropertyValue("OutputWorkspace");
int nq = static_cast<int>( qvmod.size() ); //number of q-vectors

DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", 2, nq, nq));
const std::string wsName = gwsName + std::string("_") + setName;
ws->setTitle(wsName);

double* buf = new double[nq*2];
this->dataSetDouble(h5file,setName,buf);
MantidVec& re = ws->dataY(0); // store the real part
ws->dataX(0) = qvmod; //X-axis values are the modulus of the q vector
MantidVec& im = ws->dataY(1); // store the imaginary part
ws->dataX(1) = qvmod;
double* curr = buf;
for(int iq=0; iq<nq; iq++){
re[iq]=curr[0];
im[iq]=curr[1];
curr += 2;
}
delete buf;
this->registerWorkspace(gws,wsName,ws, "X-axis: Q-vector modulus; Y-axis: intermediate structure factor");
}

/**
* Create one workspace to hold the real part and another to hold the imaginary part
* Y-values are structure factor for each Q-value
* X-values are time bins
* @param h5file file identifier
* @param gws pointer to WorkspaceGroup being filled
* @param setName string name of dataset
* @param qvmod vector of Q-vectors' moduli
*/
void LoadSassena::loadFQT(const hid_t& h5file, API::WorkspaceGroup_sptr gws, const std::string setName, const MantidVec &qvmod)
{
const std::string gwsName = this->getPropertyValue("OutputWorkspace");
int nq = static_cast<int>( qvmod.size() ); //number of q-vectors
double dt = 1.0; //time unit increment;
hsize_t dims[3];
this->dataSetInfo(h5file, setName, dims);
int nt = static_cast<int>( dims[1] ); //number of time points
double* buf = new double[nq*nt*2];
this->dataSetDouble(h5file,setName,buf);

DataObjects::Workspace2D_sptr wsRe = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
const std::string wsReName = gwsName + std::string("_") + setName + std::string(".Re");
wsRe->setTitle(wsReName);

DataObjects::Workspace2D_sptr wsIm = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
const std::string wsImName = gwsName + std::string("_") + setName + std::string(".Im");
wsIm->setTitle(wsImName);

double* curr = buf;
for(int iq=0; iq<nq; iq++)
{
MantidVec& reX = wsRe->dataX(iq);
MantidVec& imX = wsIm->dataX(iq);
MantidVec& reY = wsRe->dataY(iq);
MantidVec& imY = wsIm->dataY(iq);
for(int it=0; it<nt; it++)
{
reX[it] = it*dt; // time point for the real part
reY[it] = *curr; // real part of the intermediate structure factor
curr ++;
imX[it] = it*dt;
imY[it] = *curr;
curr ++;
}
}
delete buf;
this->registerWorkspace(gws,wsReName,wsRe, "X-axis: time; Y-axis: real part of intermediate structure factor");
this->registerWorkspace(gws,wsImName,wsIm, "X-axis: time; Y-axis: imaginary part of intermediate structure factor");
}

/**
* Initialise the algorithm. Declare properties which can be set before execution (input) or
* read from after the execution (output).
Expand Down Expand Up @@ -97,120 +260,40 @@ void LoadSassena::exec()
const char* validSets[] = { "fq", "fq0", "fq2", "fqt"};
for(int iSet=0; iSet<nvalidSets; iSet++) this->m_validSets.push_back( validSets[iSet] );

//open the HDF5 file
//open the HDF5 file for reading
m_filename = this->getPropertyValue("Filename");
hid_t h5file = H5Fopen(m_filename.c_str(),H5F_ACC_RDONLY,H5P_DEFAULT);
herr_t status;
if( h5file < 0)
{
this->g_log.error("Cannot open "+m_filename);
throw Kernel::Exception::FileError("Unable to open:" , m_filename);
}

//find out the sassena version used
char cversion[16];
status = H5LTget_attribute_string( h5file, "/", "sassena_version", cversion );
if ( H5LTget_attribute_string( h5file, "/", "sassena_version", cversion ) < 0 )
{
this->g_log.error("Cannot read version string");
throw Kernel::Exception::FileError("Unable to read version string:" , m_filename);
}
const std::string version(cversion);

//determine which loader protocol to use based on the version
//to be done at a later time, maybe implement a Version class

int rank[3]; //store dimensions
H5T_class_t class_id;
size_t type_size;
hsize_t dims[3];
/*load vectors onto a Workspace2D with 3 bins (the three components of the vectors)
* dataX for the origin of the vector (assumed (0,0,0) )
* dataY for the tip of the vector
* dataE is assumed (0,0,0), no errors
*/
std::string setName("qvectors");
status = H5LTget_dataset_info( h5file, setName.c_str(), dims, &class_id, &type_size );
int nq = static_cast<int>( dims[0] ); //number of q-vectors
double* buf = new double[nq*3];
H5LTread_dataset_double(h5file,"qvectors",buf);

DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, 3, 3));
std::string wsName = gwsName + std::string(".") + setName;
ws->setTitle(wsName);
MantidVec qvmod; //store the modulus of the vector

double* curr = buf;
for(int iq=0; iq<nq; iq++){
MantidVec& Y = ws->dataY(iq);
Y.assign(curr,curr+3);
qvmod.push_back( sqrt( curr[0]*curr[0] + curr[1]*curr[1] + curr[2]*curr[2] ) );
curr += 3;
}

delete buf;
this->registerWorkspace(gws,setName,ws, "X-axis: origin of Q-vectors; Y-axis: tip of Q-vectors");


const MantidVec qvmod = this->loadQvectors( h5file, gws );
//iterate over the valid sets
std::string setName;
for(std::vector<std::string>::const_iterator it=this->m_validSets.begin(); it!=this->m_validSets.end(); ++it){
setName = *it;
if(H5LTfind_dataset(h5file,setName.c_str())==1)
{
if(setName == "fq" || setName == "fq0" || setName == "fq2")
/* create workspace to store the structure factor.
* First spectrum is the real part, second spectrum is the imaginary part
*/
{
DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", 2, nq, nq));
wsName = gwsName + std::string(".") + setName;
ws->setTitle(wsName);
double* buf = new double[nq*2];
status = H5LTread_dataset_double(h5file,setName.c_str(),buf);
MantidVec& re = ws->dataY(0); // store the real part
ws->dataX(0) = qvmod; //X-axis values are the modulus of the q vector
MantidVec& im = ws->dataY(1); // store the imaginary part
ws->dataX(1) = qvmod;
double* curr = buf;
for(int iq=0; iq<nq; iq++){
re[iq]=curr[0];
im[iq]=curr[1];
curr += 2;
}
delete buf;
this->registerWorkspace(gws,setName,ws, "X-axis: Q-vector modulus; Y-axis: intermediate structure factor");
}

this->loadFQ( h5file, gws, setName, qvmod );
else if(setName == "fqt")
//Create one workspace to hold the real part and another to hold the imaginary part
{
double dt = 1.0; //time unit increment;
status = H5LTget_dataset_ndims( h5file, setName.c_str(), rank );
status = H5LTget_dataset_info( h5file, setName.c_str(), dims, &class_id, &type_size ); //dims = {nq,nt,2}
int nt = static_cast<int>( dims[1] ); //number of time points
double* buf = new double[nq*nt*2];
status = H5LTread_dataset_double(h5file,setName.c_str(),buf);

DataObjects::Workspace2D_sptr wsRe = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
const std::string wsReName = gwsName + std::string(".") + setName + std::string(".Re");
wsRe->setTitle(wsReName);

DataObjects::Workspace2D_sptr wsIm = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
const std::string wsImName = gwsName + std::string(".") + setName + std::string(".Im");
wsIm->setTitle(wsImName);

double* curr = buf;
for(int iq=0; iq<nq; iq++)
{
MantidVec& reX = wsRe->dataX(iq);
MantidVec& imX = wsIm->dataX(iq);
MantidVec& reY = wsRe->dataY(iq);
MantidVec& imY = wsIm->dataY(iq);
for(int it=0; it<nt; it++)
{
reX[it] = it*dt; // time point for the real part
reY[it] = *curr; // real part of the intermediate structure factor
curr ++;
imX[it] = it*dt;
imY[it] = *curr;
curr ++;
}
}
delete buf;
this->registerWorkspace(gws,wsReName,wsRe, "X-axis: time; Y-axis: real part of intermediate structure factor");
this->registerWorkspace(gws,wsImName,wsIm, "X-axis: time; Y-axis: imaginary part of intermediate structure factor");
} // end of else if(setName == "fqt")
} // end of if(H5LTfind_dataset...
this->loadFQT( h5file, gws, setName, qvmod );
}
else
this->g_log.information("Dataset "+setName+" not present in file");
}// end of iterate over the valid sets

this->setProperty( "OutputWorkspace", gws ); //register the groupWorkspace in the analysis data service
Expand Down

0 comments on commit 639561c

Please sign in to comment.