Newer
Older
Sofia Antony
committed
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/LoadRawHelper.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidAPI/MemoryManager.h"
#include "MantidAPI/SpectraDetectorMap.h"
#include "MantidAPI/WorkspaceGroup.h"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Instrument/XMLlogfile.h"
Sofia Antony
committed
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/ArrayProperty.h"
Gigg, Martyn Anthony
committed
#include "MantidAPI/FileProperty.h"
Sofia Antony
committed
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/ListValidator.h"
Sofia Antony
committed
#include "LoadRaw/isisraw2.h"
#include "MantidDataHandling/LoadLog.h"
Sofia Antony
committed
#include "MantidAPI/LoadAlgorithmFactory.h"
Robert Whitley
committed
#include "MantidDataHandling/LoadAscii.h"
Sofia Antony
committed
#include <boost/shared_ptr.hpp>
Sofia Antony
committed
#include <cmath>
#include <cstdio> //Required for gcc 4.4
#include "MantidKernel/Strings.h"
Sofia Antony
committed
namespace Mantid
{
Gigg, Martyn Anthony
committed
namespace DataHandling
{
Sofia Antony
committed
Gigg, Martyn Anthony
committed
using namespace Kernel;
using namespace API;
Sofia Antony
committed
Gigg, Martyn Anthony
committed
LoadRawHelper::LoadRawHelper() :
isisRaw(new ISISRAW2),
m_list(false),m_spec_list(),m_spec_min(0),
m_spec_max(EMPTY_INT()), m_numberOfPeriods(0), m_specTimeRegimes(),m_prog(0),m_bmspeclist(false)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
LoadRawHelper::~LoadRawHelper()
{
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/// Initialisation method.
void LoadRawHelper::init()
{
std::vector<std::string> exts;
exts.push_back(".raw");
exts.push_back(".s*");
exts.push_back(".add");
declareProperty(new FileProperty("Filename", "", FileProperty::Load, exts),
"The name of the RAW file to read, including its full or relative\n"
"path. (N.B. case sensitive if running on Linux).");
Gigg, Martyn Anthony
committed
declareProperty(new WorkspaceProperty<Workspace> ("OutputWorkspace", "", Direction::Output),
"The name of the workspace that will be created, filled with the\n"
"read-in data and stored in the Analysis Data Service. If the input\n"
"RAW file contains multiple periods higher periods will be stored in\n"
"separate workspaces called OutputWorkspace_PeriodNo.");
Gigg, Martyn Anthony
committed
Sofia Antony
committed
Gigg, Martyn Anthony
committed
m_cache_options.push_back("If Slow");
m_cache_options.push_back("Always");
m_cache_options.push_back("Never");
declareProperty("Cache", "If Slow", boost::make_shared<StringListValidator>(m_cache_options));
Sofia Antony
committed
Gigg, Martyn Anthony
committed
declareProperty("LoadLogFiles", true, " Boolean option to load or skip log files.");
Gigg, Martyn Anthony
committed
Gigg, Martyn Anthony
committed
}
/**opens the raw file and returns the file pointer
*@param fileName :: name of the raw file
*@return file pointer
*/
Gigg, Martyn Anthony
committed
FILE* LoadRawHelper::openRawFile(const std::string & fileName)
{
FILE* file = fopen(fileName.c_str(), "rb");
if (file == NULL)
{
Gigg, Martyn Anthony
committed
g_log.error("Unable to open file " + fileName);
throw Exception::FileError("Unable to open File:", fileName);
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
// Need to check that the file is not a text file as the ISISRAW routines don't deal with these very well, i.e
// reading continues until a bad_alloc is encountered.
if( isAscii(file) )
{
Gigg, Martyn Anthony
committed
g_log.error() << "File \"" << fileName << "\" is not a valid RAW file.\n";
fclose(file);
throw std::invalid_argument("Incorrect file type encountered.");
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
Gigg, Martyn Anthony
committed
return file;
Sofia Antony
committed
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
/** Reads the run title and creates a string from it
* @param file :: pointer to the raw file
* @param title :: An output parameter that will contain the workspace title
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::readTitle(FILE* file,std::string & title)
{
Gigg, Martyn Anthony
committed
ioRaw(file, true);
Gigg, Martyn Anthony
committed
title = std::string(isisRaw->r_title, 80);
Gigg, Martyn Anthony
committed
g_log.information("*** Run title: " + title + " ***");
}
/**skips the histogram from raw file
*@param file :: pointer to the raw file
*@param hist :: postion in the file to skip
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::skipData(FILE* file,int hist)
{
isisRaw->skipData(file, hist);
}
void LoadRawHelper::skipData(FILE* file,int64_t hist)
{
skipData(file, static_cast<int>(hist));
Gigg, Martyn Anthony
committed
/// calls isisRaw ioRaw.
Janik Zikovsky
committed
/// @param file :: the file pointer
/// @param from_file :: unknown
Gigg, Martyn Anthony
committed
void LoadRawHelper::ioRaw(FILE* file,bool from_file )
{
isisRaw->ioRAW(file, from_file);
}
int LoadRawHelper::getNumberofTimeRegimes()
{
return isisRaw->daep.n_tr_shift;
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
void LoadRawHelper::reset()
{
isisRaw.reset();
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/**reads the histogram from raw file
Janik Zikovsky
committed
* @param file :: pointer to the raw file
* @param hist :: postion in the file to read
* @return flag is data is read
*/
Gigg, Martyn Anthony
committed
bool LoadRawHelper::readData(FILE* file,int hist)
{
return isisRaw->readData(file, hist);
}
bool LoadRawHelper::readData(FILE* file,int64_t hist)
{
return readData(file, static_cast<int>(hist));
Sofia Antony
committed
Gigg, Martyn Anthony
committed
float LoadRawHelper::getProtonCharge()const
{
return isisRaw->rpb.r_gd_prtn_chrg;
}
Gigg, Martyn Anthony
committed
Gigg, Martyn Anthony
committed
/**
* Set the proton charge on the run object
* @param run :: The run object
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::setProtonCharge(API::Run& run)
{
run.setProtonCharge(getProtonCharge());
}
/** Stores the run number in the run logs
* @param run :: the workspace's run object
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::setRunNumber(API::Run& run)
{
std::string run_num = boost::lexical_cast<std::string>(isisRaw->r_number);
run.addLogData( new PropertyWithValue<std::string>("run_number", run_num) );
}
/**reads workspace dimensions,number of periods etc from raw data
Janik Zikovsky
committed
* @param numberOfSpectra :: number of spectrums
* @param numberOfPeriods :: number of periods
* @param lengthIn :: size of workspace vectors
* @param noTimeRegimes :: number of time regime.
void LoadRawHelper::readworkspaceParameters(specid_t& numberOfSpectra,int& numberOfPeriods,int64_t& lengthIn,int64_t & noTimeRegimes )
Gigg, Martyn Anthony
committed
{
// Read in the number of spectra in the RAW file
m_numberOfSpectra=numberOfSpectra = static_cast<specid_t>(isisRaw->t_nsp1);
Gigg, Martyn Anthony
committed
// Read the number of periods in this file
numberOfPeriods = isisRaw->t_nper;
// Read the number of time channels (i.e. bins) from the RAW file
Gigg, Martyn Anthony
committed
// Read in the time bin boundaries
lengthIn = channelsPerSpectrum + 1;
// Now check whether there is more than one time regime in use
noTimeRegimes = isisRaw->daep.n_tr_shift;
}
/**This method creates shared pointer to a workspace
Janik Zikovsky
committed
* @param ws_sptr :: shared pointer to the parent workspace
* @param nVectors :: number of histograms in the workspace
* @param xLengthIn :: size of workspace X vector
* @param yLengthIn :: size of workspace Y vector
* @return an empty workspace of the given parameters
*/
Gigg, Martyn Anthony
committed
DataObjects::Workspace2D_sptr LoadRawHelper::createWorkspace(DataObjects::Workspace2D_sptr ws_sptr,
int64_t nVectors,int64_t xLengthIn,int64_t yLengthIn)
Gigg, Martyn Anthony
committed
{
DataObjects::Workspace2D_sptr empty;
if(!ws_sptr)return empty;
DataObjects::Workspace2D_sptr workspace = boost::dynamic_pointer_cast<DataObjects::Workspace2D>
(WorkspaceFactory::Instance().create(ws_sptr,nVectors,xLengthIn,yLengthIn));
Gigg, Martyn Anthony
committed
return workspace;
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** This method creates pointer to workspace
* @param nVectors :: The number of vectors/histograms in the workspace
* @param xlengthIn :: The number of X data points/bin boundaries in each vector
* @param ylengthIn :: The number of Y data points/bin boundaries in each vector
* @param title :: title of the workspace
* @return Workspace2D_sptr shared pointer to the workspace
*/
DataObjects::Workspace2D_sptr LoadRawHelper::createWorkspace(int64_t nVectors, int64_t xlengthIn,int64_t ylengthIn,const std::string& title)
Gigg, Martyn Anthony
committed
{
DataObjects::Workspace2D_sptr workspace;
if(nVectors>0)
{
Gigg, Martyn Anthony
committed
workspace =boost::dynamic_pointer_cast<DataObjects::Workspace2D>(WorkspaceFactory::Instance().create(
"Workspace2D", nVectors, xlengthIn, ylengthIn));
Gigg, Martyn Anthony
committed
// Set the units
workspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
workspace->setYUnit("Counts");
workspace->setTitle(title);
Gigg, Martyn Anthony
committed
}
return workspace;
}
/**creates monitor workspace
*@param monws_sptr :: shared pointer to monitor workspace
*@param normalws_sptr :: shared pointer to output workspace
*@param mongrp_sptr :: shared pointer to monitor group workspace
*@param mwsSpecs :: number of spectra in the monitor workspace
*@param nwsSpecs :: number of spectra in the output workspace
*@param numberOfPeriods :: total number of periods from raw file
*@param lengthIn :: size of workspace vectors
*@param title :: title of the workspace
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::createMonitorWorkspace(DataObjects::Workspace2D_sptr& monws_sptr,DataObjects::Workspace2D_sptr& normalws_sptr,
WorkspaceGroup_sptr& mongrp_sptr,const int64_t mwsSpecs,const int64_t nwsSpecs,
const int64_t numberOfPeriods,const int64_t lengthIn,const std::string title)
Gigg, Martyn Anthony
committed
{
try
{
Gigg, Martyn Anthony
committed
//create monitor group workspace
mongrp_sptr = createGroupWorkspace(); //create workspace
// create monitor workspace
if(mwsSpecs>0)
{
if(normalws_sptr)
{
monws_sptr=createWorkspace(normalws_sptr,mwsSpecs,lengthIn,lengthIn-1);
}
else
{
monws_sptr=createWorkspace(mwsSpecs,lengthIn,lengthIn-1,title);
}
}
if(!monws_sptr) return ;
std::string wsName= getPropertyValue("OutputWorkspace");
// if the normal output workspace size>0 then set the workspace as "MonitorWorkspace"
// otherwise set the workspace as "OutputWorkspace"
if (nwsSpecs> 0)
{
std::string monitorwsName = wsName + "_Monitors";
declareProperty(new WorkspaceProperty<Workspace> ("MonitorWorkspace", monitorwsName,
Gigg, Martyn Anthony
committed
setWorkspaceProperty("MonitorWorkspace", title, mongrp_sptr, monws_sptr,numberOfPeriods, true);
}
else
{
//if only monitors range selected
//then set the monitor workspace as the outputworkspace
setWorkspaceProperty("OutputWorkspace", title, mongrp_sptr, monws_sptr,numberOfPeriods, false);
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
}
catch(std::out_of_range& )
{
Gigg, Martyn Anthony
committed
g_log.debug()<<"Error in creating monitor workspace"<<std::endl;
Gigg, Martyn Anthony
committed
}
catch(std::runtime_error& )
{
Gigg, Martyn Anthony
committed
g_log.debug()<<"Error in creating monitor workspace"<<std::endl;
Gigg, Martyn Anthony
committed
}
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** Executes the algorithm. Reading in the file and creating and populating
* the output workspace
*
* @throw Exception::FileError If the RAW file cannot be found/opened
* @throw std::invalid_argument If the optional properties are set to invalid values
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::exec()
Sofia Antony
committed
{
}
Gigg, Martyn Anthony
committed
/** sets the workspace properties
* @param ws_sptr :: shared pointer to workspace
* @param grpws_sptr :: shared pointer to group workspace
* @param period period number
* @param bmonitors :: boolean flag to name the workspaces
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::setWorkspaceProperty(DataObjects::Workspace2D_sptr ws_sptr, WorkspaceGroup_sptr grpws_sptr,
const int64_t period, bool bmonitors)
Gigg, Martyn Anthony
committed
{
if(!ws_sptr) return;
if(!grpws_sptr) return;
std::string wsName;
std::string outws;
std::string outputWorkspace;
std::string localWSName = getProperty("OutputWorkspace");
std::stringstream suffix;
suffix << (period + 1);
if (bmonitors)
{
Gigg, Martyn Anthony
committed
wsName = localWSName + "_Monitors" + "_" + suffix.str();
outputWorkspace = "MonitorWorkspace";
Gigg, Martyn Anthony
committed
}
else
{
Gigg, Martyn Anthony
committed
wsName = localWSName + "_" + suffix.str();
outputWorkspace = "OutputWorkspace";
Gigg, Martyn Anthony
committed
}
outws = outputWorkspace + "_" + suffix.str();
declareProperty(new WorkspaceProperty<Workspace> (outws, wsName, Direction::Output));
setProperty(outws, boost::static_pointer_cast<Workspace>(ws_sptr));
grpws_sptr->addWorkspace( ws_sptr );
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** This method sets the workspace property
* @param propertyName :: property name for the workspace
* @param title :: title of the workspace
* @param grpws_sptr :: shared pointer to group workspace
* @param ws_sptr :: shared pointer to workspace
* @param numberOfPeriods :: numer periods in the raw file
* @param bMonitor to identify the workspace is an output workspace or monitor workspace
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::setWorkspaceProperty(const std::string& propertyName, const std::string& title,
WorkspaceGroup_sptr grpws_sptr, DataObjects::Workspace2D_sptr ws_sptr,int64_t numberOfPeriods, bool bMonitor)
Gigg, Martyn Anthony
committed
{
Gigg, Martyn Anthony
committed
Property *ws = getProperty("OutputWorkspace");
if(!ws) return;
if(!grpws_sptr) return;
if(!ws_sptr)return;
Gigg, Martyn Anthony
committed
ws_sptr->setTitle(title);
ws_sptr->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
if (numberOfPeriods > 1)
{
Gigg, Martyn Anthony
committed
setProperty(propertyName, boost::dynamic_pointer_cast<Workspace>(grpws_sptr));
Gigg, Martyn Anthony
committed
}
else
{
Gigg, Martyn Anthony
committed
setProperty(propertyName, boost::dynamic_pointer_cast<Workspace>(ws_sptr));
Gigg, Martyn Anthony
committed
}
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** This method sets the raw file data to workspace vectors
* @param newWorkspace :: shared pointer to the workspace
* @param timeChannelsVec :: vector holding the X data
* @param wsIndex variable used for indexing the ouputworkspace
* @param nspecNum spectrum number
* @param noTimeRegimes :: regime no.
* @param lengthIn :: length of the workspace
* @param binStart :: start of bin
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::setWorkspaceData(DataObjects::Workspace2D_sptr newWorkspace, const std::vector<
boost::shared_ptr<MantidVec> >& timeChannelsVec, int64_t wsIndex, specid_t nspecNum, int64_t noTimeRegimes,int64_t lengthIn,int64_t binStart)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
if(!newWorkspace)return;
typedef double (*uf)(double);
uf dblSqrt = std::sqrt;
// But note that the last (overflow) bin is kept
MantidVec& Y = newWorkspace->dataY(wsIndex);
Y.assign(isisRaw->dat1 + binStart, isisRaw->dat1 + lengthIn);
// Fill the vector for the errors, containing sqrt(count)
MantidVec& E = newWorkspace->dataE(wsIndex);
std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
newWorkspace->getAxis(1)->spectraNo(wsIndex) = nspecNum;
//for loadrawbin0
if(binStart==0)
{
Gigg, Martyn Anthony
committed
newWorkspace->setX(wsIndex, timeChannelsVec[0]);
return ;
Gigg, Martyn Anthony
committed
}
//for loadrawspectrum 0
if(nspecNum==0)
{
Gigg, Martyn Anthony
committed
newWorkspace->setX(wsIndex, timeChannelsVec[0]);
return ;
Gigg, Martyn Anthony
committed
}
// Set the X vector pointer and spectrum number
if (noTimeRegimes < 2)
Gigg, Martyn Anthony
committed
newWorkspace->setX(wsIndex, timeChannelsVec[0]);
Gigg, Martyn Anthony
committed
else
{
Gigg, Martyn Anthony
committed
// Use std::vector::at just incase spectrum missing from spec array
newWorkspace->setX(wsIndex, timeChannelsVec.at(m_specTimeRegimes[nspecNum] - 1));
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
/** This method returns the monitor spectrum list
* @param localWorkspace :: shared pointer to workspace
* @param monitorSpecList :: a list holding the spectrum indexes of the monitors
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::getmonitorSpectrumList(DataObjects::Workspace2D_sptr localWorkspace,
std::vector<specid_t>& monitorSpecList)
Gigg, Martyn Anthony
committed
{
if (!m_monitordetectorList.empty())
{
std::vector<specid_t> specList;
//get the monitor spectrum list from SpectraDetectorMap
localWorkspace->getSpectraFromDetectorIDs(m_monitordetectorList, specList);
// Old way to get the spectra # for these detectors
Gigg, Martyn Anthony
committed
const Geometry::ISpectraDetectorMap& specdetMap = localWorkspace->spectraMap();
specList = specdetMap.getSpectra(m_monitordetectorList);
Gigg, Martyn Anthony
committed
// remove duplicates by calling sort & unique algorithms
sort(specList.begin(), specList.end(), std::less<int>());
std::vector<specid_t>::iterator uEnd;
Gigg, Martyn Anthony
committed
uEnd = unique(specList.begin(), specList.end());
std::vector<specid_t> newVec;
Gigg, Martyn Anthony
committed
newVec.assign(specList.begin(), uEnd);
//remove if zeroes are there in the Spectra list
std::vector<specid_t>::iterator itr;
Gigg, Martyn Anthony
committed
itr = find(newVec.begin(), newVec.end(), 0);
if (itr != newVec.end())
newVec.erase(itr);
monitorSpecList.assign(newVec.begin(), newVec.end());
Gigg, Martyn Anthony
committed
}
else{
Gigg, Martyn Anthony
committed
g_log.error() << "monitor detector id list is empty for the selected workspace" << std::endl;
Gigg, Martyn Anthony
committed
}
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** This method creates pointer to group workspace
* @return WorkspaceGroup_sptr shared pointer to the workspace
*/
Gigg, Martyn Anthony
committed
WorkspaceGroup_sptr LoadRawHelper::createGroupWorkspace()
{
WorkspaceGroup_sptr workspacegrp(new WorkspaceGroup);
return workspacegrp;
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/**
* Check if a file is a text file
* @param file :: The file pointer
* @returns true if the file an ascii text file, false otherwise
*/
Gigg, Martyn Anthony
committed
bool LoadRawHelper::isAscii(FILE* file) const
{
Robert Whitley
committed
return LoadAscii::isAscii(file);
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
Gigg, Martyn Anthony
committed
/** Constructs the time channel (X) vector(s)
* @param regimes :: The number of time regimes (if 1 regime, will actually contain 0)
* @param lengthIn :: The number of time channels
* @return The vector(s) containing the time channel boundaries, in a vector of shared ptrs
*/
std::vector<boost::shared_ptr<MantidVec> > LoadRawHelper::getTimeChannels(const int64_t& regimes,
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
float* const timeChannels = new float[lengthIn];
isisRaw->getTimeChannels(timeChannels, static_cast<int>(lengthIn));
Gigg, Martyn Anthony
committed
std::vector<boost::shared_ptr<MantidVec> > timeChannelsVec;
if (regimes >= 2)
{
Gigg, Martyn Anthony
committed
g_log.debug() << "Raw file contains " << regimes << " time regimes\n";
// If more than 1 regime, create a timeChannelsVec for each regime
Gigg, Martyn Anthony
committed
{
// Create a vector with the 'base' time channels
boost::shared_ptr<MantidVec> channelsVec(new MantidVec(timeChannels, timeChannels + lengthIn));
const double shift = isisRaw->daep.tr_shift[i];
g_log.debug() << "Time regime " << i + 1 << " shifted by " << shift << " microseconds\n";
// Add on the shift for this vector
std::transform(channelsVec->begin(), channelsVec->end(), channelsVec->begin(), std::bind2nd(
std::plus<double>(), shift));
Gigg, Martyn Anthony
committed
timeChannelsVec.push_back(channelsVec);
}
// In this case, also need to populate the map of spectrum-regime correspondence
const int64_t ndet = static_cast<int64_t>(isisRaw->i_det);
Janik Zikovsky
committed
std::map<specid_t, specid_t>::iterator hint = m_specTimeRegimes.begin();
Gigg, Martyn Anthony
committed
{
// No checking for consistency here - that all detectors for given spectrum
// are declared to use same time regime. Will just use first encountered
hint = m_specTimeRegimes.insert(hint, std::make_pair(isisRaw->spec[j], isisRaw->timr[j]));
}
Gigg, Martyn Anthony
committed
}
else // Just need one in this case
{
Gigg, Martyn Anthony
committed
boost::shared_ptr<MantidVec> channelsVec(new MantidVec(timeChannels, timeChannels + lengthIn));
timeChannelsVec.push_back(channelsVec);
Gigg, Martyn Anthony
committed
}
// Done with the timeChannels C array so clean up
delete[] timeChannels;
return timeChannelsVec;
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
/// Run the sub-algorithm LoadInstrument (or LoadInstrumentFromRaw)
Janik Zikovsky
committed
/// @param fileName :: the raw file filename
/// @param localWorkspace :: The workspace to load the instrument for
/// @param progStart :: progress at start
/// @param progEnd :: progress at end
void LoadRawHelper::runLoadInstrument(const std::string& fileName,DataObjects::Workspace2D_sptr localWorkspace, double progStart, double progEnd)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
g_log.debug("Loading the instrument definition...");
m_prog = progStart;
Gigg, Martyn Anthony
committed
progress(m_prog, "Loading the instrument geometry...");
std::string instrumentID = isisRaw->i_inst; // get the instrument name
size_t i = instrumentID.find_first_of(' '); // cut trailing spaces
Gigg, Martyn Anthony
committed
if (i != std::string::npos) instrumentID.erase(i);
Gigg, Martyn Anthony
committed
IAlgorithm_sptr loadInst= createSubAlgorithm("LoadInstrument");
// Enable progress reporting by sub-algorithm -
loadInst->addObserver(m_progressObserver);
setChildStartProgress(progStart);
setChildEndProgress((progStart+progEnd)/2);
Gigg, Martyn Anthony
committed
// Now execute the sub-algorithm. Catch and log any error, but don't stop.
bool executionSuccessful(true);
try
{
loadInst->setPropertyValue("InstrumentName", instrumentID);
Gigg, Martyn Anthony
committed
loadInst->setProperty<MatrixWorkspace_sptr> ("Workspace", localWorkspace);
Gigg, Martyn Anthony
committed
loadInst->setProperty("RewriteSpectraMap", false); // No point as we will load the one from the file
Gigg, Martyn Anthony
committed
loadInst->execute();
Gigg, Martyn Anthony
committed
} catch (std::invalid_argument&)
{
Gigg, Martyn Anthony
committed
g_log.information("Invalid argument to LoadInstrument sub-algorithm");
executionSuccessful = false;
Gigg, Martyn Anthony
committed
} catch (std::runtime_error&)
{
Gigg, Martyn Anthony
committed
g_log.information("Unable to successfully run LoadInstrument sub-algorithm");
executionSuccessful = false;
Gigg, Martyn Anthony
committed
}
// If loading instrument definition file fails, run LoadInstrumentFromRaw instead
if (!executionSuccessful)
{
Gigg, Martyn Anthony
committed
g_log.information() << "Instrument definition file "
<< " not found. Attempt to load information about \n"
<< "the instrument from raw data file.\n";
Gigg, Martyn Anthony
committed
runLoadInstrumentFromRaw(fileName,localWorkspace);
Gigg, Martyn Anthony
committed
}
else
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
// If requested update the instrument to positions in the raw file
const Geometry::ParameterMap & pmap = localWorkspace->instrumentParameters();
Russell Taylor
committed
if( pmap.contains(localWorkspace->getInstrument()->getComponentID(),"det-pos-source") )
Gigg, Martyn Anthony
committed
{
Russell Taylor
committed
boost::shared_ptr<Geometry::Parameter> updateDets = pmap.get(localWorkspace->getInstrument()->getComponentID(),"det-pos-source");
Roman Tolchenov
committed
std::string value = updateDets->value<std::string>();
if(value.substr(0,8) == "datafile" )
Gigg, Martyn Anthony
committed
{
IAlgorithm_sptr updateInst = createSubAlgorithm("UpdateInstrumentFromFile");
updateInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace);
updateInst->setPropertyValue("Filename", fileName);
updateInst->addObserver(m_progressObserver); // Enable progress reporting by subalgorithm
setChildStartProgress((progStart+progEnd)/2);
setChildEndProgress(progEnd);
Roman Tolchenov
committed
if(value == "datafile-ignore-phi" )
{
updateInst->setProperty("IgnorePhi", true);
g_log.information("Detector positions in IDF updated with positions in the data file except for the phi values");
}
else
{
g_log.information("Detector positions in IDF updated with positions in the data file");
Roman Tolchenov
committed
}
Gigg, Martyn Anthony
committed
// We want this to throw if it fails to warn the user that the information is not correct.
updateInst->execute();
}
}
// Debugging code??
Gigg, Martyn Anthony
committed
m_monitordetectorList = loadInst->getProperty("MonitorList");
Gigg, Martyn Anthony
committed
for (itr = m_monitordetectorList.begin(); itr != m_monitordetectorList.end(); ++itr)
{
g_log.debug() << "Monitor detector id is " << (*itr) << std::endl;
}
Sofia Antony
committed
}
}
Gigg, Martyn Anthony
committed
/// Run LoadInstrumentFromRaw as a sub-algorithm (only if loading from instrument definition file fails)
Janik Zikovsky
committed
/// @param fileName :: the raw file filename
/// @param localWorkspace :: The workspace to load the instrument for
Gigg, Martyn Anthony
committed
void LoadRawHelper::runLoadInstrumentFromRaw(const std::string& fileName,DataObjects::Workspace2D_sptr localWorkspace)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
IAlgorithm_sptr loadInst = createSubAlgorithm("LoadInstrumentFromRaw");
loadInst->setPropertyValue("Filename", fileName);
// Set the workspace property to be the same one filled above
loadInst->setProperty<MatrixWorkspace_sptr> ("Workspace", localWorkspace);
// Now execute the sub-algorithm. Catch and log any error, but don't stop.
try
{
Gigg, Martyn Anthony
committed
loadInst->execute();
Gigg, Martyn Anthony
committed
} catch (std::runtime_error&)
{
Gigg, Martyn Anthony
committed
g_log.error("Unable to successfully run LoadInstrumentFromRaw sub-algorithm");
Gigg, Martyn Anthony
committed
}
m_monitordetectorList = loadInst->getProperty("MonitorList");
Gigg, Martyn Anthony
committed
for (itr = m_monitordetectorList.begin(); itr != m_monitordetectorList.end(); ++itr)
{
Gigg, Martyn Anthony
committed
g_log.debug() << "Monitor dtector id is " << (*itr) << std::endl;
Gigg, Martyn Anthony
committed
}
if (!loadInst->isExecuted())
{
Gigg, Martyn Anthony
committed
g_log.error("No instrument definition loaded");
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
/// Run the LoadMappingTable sub-algorithm to fill the SpectraToDetectorMap
Janik Zikovsky
committed
/// @param fileName :: the raw file filename
/// @param localWorkspace :: The workspace to load the mapping table for
Gigg, Martyn Anthony
committed
void LoadRawHelper::runLoadMappingTable(const std::string& fileName,DataObjects::Workspace2D_sptr localWorkspace)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
g_log.debug("Loading the spectra-detector mapping...");
progress(m_prog, "Loading the spectra-detector mapping...");
// Now determine the spectra to detector map calling sub-algorithm LoadMappingTable
// There is a small penalty in re-opening the raw file but nothing major.
IAlgorithm_sptr loadmap = createSubAlgorithm("LoadMappingTable");
loadmap->setPropertyValue("Filename",fileName);
loadmap->setProperty<MatrixWorkspace_sptr> ("Workspace", localWorkspace);
try
{
Gigg, Martyn Anthony
committed
loadmap->execute();
Gigg, Martyn Anthony
committed
} catch (std::runtime_error&)
{
Gigg, Martyn Anthony
committed
g_log.error("Unable to successfully execute LoadMappingTable sub-algorithm");
Gigg, Martyn Anthony
committed
}
if (!loadmap->isExecuted())
{
Gigg, Martyn Anthony
committed
g_log.error("LoadMappingTable sub-algorithm is not executed");
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
/// Run the LoadLog sub-algorithm
Janik Zikovsky
committed
/// @param fileName :: the raw file filename
/// @param localWorkspace :: The workspace to load the logs for
void LoadRawHelper::runLoadLog(const std::string& fileName, DataObjects::Workspace2D_sptr localWorkspace, double progStart, double progEnd )
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
g_log.debug("Loading the log files...");
if( progStart < progEnd ) {
Gigg, Martyn Anthony
committed
progress(m_prog, "Reading log files...");
IAlgorithm_sptr loadLog = createSubAlgorithm("LoadLog");
// Pass through the same input filename
loadLog->setPropertyValue("Filename", fileName);
// Set the workspace property to be the same one filled above
loadLog->setProperty<MatrixWorkspace_sptr> ("Workspace", localWorkspace);
// Enable progress reporting by sub-algorithm - if progress range has duration
if( progStart < progEnd ) {
loadLog->addObserver(m_progressObserver);
setChildStartProgress(progStart);
setChildEndProgress(progEnd);
}
Gigg, Martyn Anthony
committed
// Now execute the sub-algorithm. Catch and log any error, but don't stop.
try
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
loadLog->execute();
Gigg, Martyn Anthony
committed
} catch (std::exception&)
{
Gigg, Martyn Anthony
committed
g_log.error("Unable to successfully run LoadLog sub-algorithm");
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
if (!loadLog->isExecuted())
{
Gigg, Martyn Anthony
committed
g_log.error("Unable to successfully run LoadLog sub-algorithm");
Gigg, Martyn Anthony
committed
}
// Make log creator object and add the run status log if we have the appropriate ICP log
m_logCreator.reset(new ISISRunLogs(localWorkspace->run(), m_numberOfPeriods));
m_logCreator->addStatusLog(localWorkspace->mutableRun());
Sofia Antony
committed
}
Gigg, Martyn Anthony
committed
/**
* Creates period log data in the workspace
* @param period :: period number
* @param local_workspace :: workspace to add period log data to.
void LoadRawHelper::createPeriodLogs(int64_t period, DataObjects::Workspace2D_sptr local_workspace)
{
m_logCreator->addPeriodLogs(static_cast<int>(period), local_workspace->mutableRun());
}
Gigg, Martyn Anthony
committed
/**
* Pulls the run parameters from the ISIS Raw RPB structure and stores them as log entries on the
* workspace run object
* @param localWorkspace :: The workspace to attach the information to
* @param rawFile :: The handle to an ISIS Raw file
*/
Gigg, Martyn Anthony
committed
void LoadRawHelper::loadRunParameters(API::MatrixWorkspace_sptr localWorkspace, ISISRAW * const rawFile) const
Gigg, Martyn Anthony
committed
{
Gigg, Martyn Anthony
committed
ISISRAW * localISISRaw(NULL);
if( !rawFile )
{
Gigg, Martyn Anthony
committed
localISISRaw = isisRaw.get();
Gigg, Martyn Anthony
committed
}
else
{
Gigg, Martyn Anthony
committed
localISISRaw = rawFile;
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
API::Run& runDetails = localWorkspace->mutableRun();
Gigg, Martyn Anthony
committed
// Run header is stored as consecutive char arrays adding up to a total of 80 bytes in the HDR_STRUCT
const std::string run_header(localISISRaw->hdr.inst_abrv, 80);
Gigg, Martyn Anthony
committed
runDetails.addProperty("run_header", run_header);
// Run title is stored in a different attribute
Anders Markvardsen
committed
runDetails.addProperty("run_title", std::string(localISISRaw->r_title,80), true);
Gigg, Martyn Anthony
committed
runDetails.addProperty("user_name", std::string(localISISRaw->hdr.hd_user, 20));
runDetails.addProperty("inst_abrv", std::string(localISISRaw->hdr.inst_abrv, 3));
runDetails.addProperty("hd_dur", std::string(localISISRaw->hdr.hd_dur, 8));
Gigg, Martyn Anthony
committed
// Data details on run not the workspace
Gigg, Martyn Anthony
committed
runDetails.addProperty("nspectra", static_cast<int>(localISISRaw->t_nsp1));
runDetails.addProperty("nchannels", static_cast<int>(localISISRaw->t_ntc1));
runDetails.addProperty("nperiods", static_cast<int>(localISISRaw->t_nper));
Gigg, Martyn Anthony
committed
// RPB struct info
Gigg, Martyn Anthony
committed
runDetails.addProperty("dur", localISISRaw->rpb.r_dur); // actual run duration
runDetails.addProperty("durunits", localISISRaw->rpb.r_durunits); // scaler for above (1=seconds)
Gigg, Martyn Anthony
committed
runDetails.addProperty("dur_freq",localISISRaw->rpb.r_dur_freq); // testinterval for above (seconds)
runDetails.addProperty("dmp", localISISRaw->rpb.r_dmp); // dump interval
Gigg, Martyn Anthony
committed
runDetails.addProperty("dmp_units", localISISRaw->rpb.r_dmp_units); // scaler for above
runDetails.addProperty("dmp_freq",localISISRaw->rpb.r_dmp_freq); // interval for above
runDetails.addProperty("freq",localISISRaw->rpb.r_freq); // 2**k where source frequency = 50 / 2**k
Gigg, Martyn Anthony
committed
runDetails.addProperty("gd_prtn_chrg", static_cast<double>(localISISRaw->rpb.r_gd_prtn_chrg)); // good proton charge (uA.hour)
runDetails.addProperty("tot_prtn_chrg", static_cast<double>(localISISRaw->rpb.r_tot_prtn_chrg)); // total proton charge (uA.hour)
Gigg, Martyn Anthony
committed
runDetails.addProperty("goodfrm",localISISRaw->rpb.r_goodfrm); // good frames
runDetails.addProperty("rawfrm", localISISRaw->rpb.r_rawfrm); // raw frames
Gigg, Martyn Anthony
committed
runDetails.addProperty("dur_wanted",localISISRaw->rpb.r_dur_wanted); // requested run duration (units as for "duration" above)
Gigg, Martyn Anthony
committed
runDetails.addProperty("dur_secs",localISISRaw->rpb.r_dur_secs ); // actual run duration in seconds
runDetails.addProperty("mon_sum1", localISISRaw->rpb.r_mon_sum1); // monitor sum 1
runDetails.addProperty("mon_sum2",localISISRaw->rpb.r_mon_sum2); // monitor sum 2
runDetails.addProperty("mon_sum3",localISISRaw->rpb.r_mon_sum3); // monitor sum 3
Gigg, Martyn Anthony
committed
runDetails.addProperty("rb_proposal",localISISRaw->rpb.r_prop); // RB (proposal) number
Gigg, Martyn Anthony
committed
// Note isis raw date format which is stored in DD-MMM-YYYY. Store dates in ISO 8601
std::string isisDate = std::string(localISISRaw->rpb.r_enddate, 11);
if ( isisDate[0] == ' ' ) isisDate[0] = '0';
runDetails.addProperty("run_end", DateAndTime(isisDate.substr(7,4) + "-" + convertMonthLabelToIntStr(isisDate.substr(3,3))
+ "-" + isisDate.substr(0,2) + "T" + std::string(localISISRaw->rpb.r_endtime, 8)).toISO8601String());
isisDate = std::string(localISISRaw->hdr.hd_date, 11);
if ( isisDate[0] == ' ' ) isisDate[0] = '0';
runDetails.addProperty("run_start", DateAndTime(isisDate.substr(7,4) + "-" + convertMonthLabelToIntStr(isisDate.substr(3,3))
+ "-" + isisDate.substr(0,2) + "T" + std::string(localISISRaw->hdr.hd_time, 8)).toISO8601String());
}
/// To help transforming date stored in ISIS raw file into iso 8601
/// @param month
:: /// @return month as string integer e.g. 01
std::string LoadRawHelper::convertMonthLabelToIntStr(std::string month) const
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
{
std::transform(month.begin(), month.end(), month.begin(), toupper);
if ( month == "JAN" )
return "01";
if ( month == "FEB" )
return "02";
if ( month == "MAR" )
return "03";
if ( month == "APR" )
return "04";
if ( month == "MAY" )
return "05";
if ( month == "JUN" )
return "06";
if ( month == "JUL" )
return "07";
if ( month == "AUG" )
return "08";
if ( month == "SEP" )
return "09";
if ( month == "OCT" )
return "10";
if ( month == "NOV" )
return "11";
if ( month == "DEC" )
return "12";
throw std::runtime_error("LoadRawHelper::convertMonthLabelToIntStr(): Invalid month label found.");
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
///sets optional properties for the loadraw algorithm
Janik Zikovsky
committed
/// @param spec_min :: The minimum spectra number
/// @param spec_max :: The maximum spectra number
/// @param spec_list :: The list of Spectra to be included
Gigg, Martyn Anthony
committed
void LoadRawHelper::setOptionalProperties(const int& spec_min,const int& spec_max,const std::vector<int>& spec_list)
{
m_spec_min=spec_min;
m_spec_max=spec_max;
m_spec_list.assign(spec_list.begin(),spec_list.end());
}
/// Validates the optional 'spectra to read' properties, if they have been set
void LoadRawHelper::checkOptionalProperties()
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
//read in the settings passed to the algorithm
/*m_spec_list = getProperty("SpectrumList");
Gigg, Martyn Anthony
committed
m_spec_max = getProperty("SpectrumMax");
m_spec_min = getProperty("SpectrumMin");*/
Gigg, Martyn Anthony
committed
m_list = !m_spec_list.empty();
m_bmspeclist = !m_spec_list.empty();
m_interval = (m_spec_max != EMPTY_INT()) || (m_spec_min != 1);
if (m_spec_max == EMPTY_INT())
Gigg, Martyn Anthony
committed
m_spec_max = 1;
Gigg, Martyn Anthony
committed
// Check validity of spectra list property, if set
if (m_list)
{
Gigg, Martyn Anthony
committed
m_list = true;
if (m_spec_list.size() == 0)
{
m_list = false;
}
else
{
const int64_t minlist = *min_element(m_spec_list.begin(), m_spec_list.end());
const int64_t maxlist = *max_element(m_spec_list.begin(), m_spec_list.end());
Gigg, Martyn Anthony
committed
if (maxlist >m_numberOfSpectra || minlist <= 0)
{
g_log.error("Invalid list of spectra");
throw std::invalid_argument("Inconsistent properties defined");
}
}
Gigg, Martyn Anthony
committed
}
// Check validity of spectra range, if set
if (m_interval)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
m_interval = true;
//m_spec_min = getProperty("SpectrumMin");
if (m_spec_min != 1 && m_spec_max == 1)
m_spec_max = m_numberOfSpectra;
if (m_spec_max < m_spec_min || m_spec_max >m_numberOfSpectra)
{
g_log.error("Invalid Spectrum min/max properties");
throw std::invalid_argument("Inconsistent properties defined");
}
Gigg, Martyn Anthony
committed
Sofia Antony
committed
}
}
/**
* Calculates the total number of spectra in the workspace, given the input properties
* @return the size of the workspace (number of spectra)
*/
specid_t LoadRawHelper::calculateWorkspaceSize()
Sofia Antony
committed
{
specid_t total_specs(0);
Gigg, Martyn Anthony
committed
if (m_interval || m_list)
Sofia Antony
committed
{
Gigg, Martyn Anthony
committed
if (m_interval)
{
if (m_spec_min != 1 && m_spec_max == 1)
m_spec_max = m_numberOfSpectra;
m_total_specs=total_specs = (m_spec_max - m_spec_min + 1);
m_spec_max += 1;
}
else
total_specs = 0;
if (m_list)
{
if (m_interval)
{
Janik Zikovsky
committed
for (std::vector<specid_t>::iterator it = m_spec_list.begin(); it != m_spec_list.end();)
Gigg, Martyn Anthony
committed
if (*it >= m_spec_min && *it < m_spec_max)
{
it = m_spec_list.erase(it);
}
else
Gigg, Martyn Anthony
committed
}
if (m_spec_list.size() == 0)
m_list = false;
total_specs += static_cast<specid_t>(m_spec_list.size());
Gigg, Martyn Anthony
committed
m_total_specs=total_specs;
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
Sofia Antony
committed
else
Gigg, Martyn Anthony
committed
{
Gigg, Martyn Anthony
committed
total_specs = m_numberOfSpectra;
Gigg, Martyn Anthony
committed
m_total_specs=total_specs;
Gigg, Martyn Anthony
committed
// In this case want all the spectra, but zeroth spectrum is garbage so go from 1 to NSP1
m_spec_min = 1;
m_spec_max = m_numberOfSpectra + 1;
Gigg, Martyn Anthony
committed
}
return total_specs;
}
/// calculate workspace sizes.
Janik Zikovsky
committed
/// @param monitorSpecList :: the vector of the monitor spectra
/// @param normalwsSpecs :: the spectra for the detector workspace
/// @param monitorwsSpecs :: the spectra for the monitor workspace
void LoadRawHelper::calculateWorkspacesizes(const std::vector<specid_t>& monitorSpecList,
specid_t& normalwsSpecs, specid_t & monitorwsSpecs)
Gigg, Martyn Anthony
committed
{
if (!m_interval && !m_bmspeclist)
{
monitorwsSpecs = static_cast<specid_t>(monitorSpecList.size());
normalwsSpecs = m_total_specs - monitorwsSpecs;
Gigg, Martyn Anthony
committed
g_log.debug() << "normalwsSpecs when m_interval & m_bmspeclist are false is " << normalwsSpecs
<< " monitorwsSpecs is " << monitorwsSpecs << std::endl;
Gigg, Martyn Anthony
committed
}
else if (m_interval || m_bmspeclist)
{
Gigg, Martyn Anthony
committed
if (m_interval)
{
Gigg, Martyn Anthony
committed
for (itr1 = monitorSpecList.begin(); itr1 != monitorSpecList.end(); ++itr1)
{
if (*itr1 >= m_spec_min && *itr1 < m_spec_max)
++msize;
}
monitorwsSpecs = msize;
normalwsSpecs = m_total_specs - monitorwsSpecs;
g_log.debug() << "normalwsSpecs when m_interval true is " << normalwsSpecs
<< " monitorwsSpecs is " << monitorwsSpecs << std::endl;
Gigg, Martyn Anthony
committed
}
if (m_bmspeclist)
{
if (m_interval)
{
Gigg, Martyn Anthony
committed
for (itr = m_spec_list.begin(); itr != m_spec_list.end();)
{ //if the m_spec_list elements are in the range between m_spec_min & m_spec_max
if (*itr >= m_spec_min && *itr < m_spec_max)
itr = m_spec_list.erase(itr);
else
++itr;
}
if (m_spec_list.size() == 0)
{
g_log.debug() << "normalwsSpecs is " << normalwsSpecs << " monitorwsSpecs is "
<< monitorwsSpecs << std::endl;
Gigg, Martyn Anthony
committed
}
else
{ //at this point there are monitors in the list which are not in the min& max range
// so find those monitors count and calculate the workspace specs
std::vector<specid_t>::const_iterator itr;
std::vector<specid_t>::const_iterator monitr;
specid_t monCounter = 0;
Gigg, Martyn Anthony
committed
for (itr = m_spec_list.begin(); itr != m_spec_list.end(); ++itr)
{
monitr = find(monitorSpecList.begin(), monitorSpecList.end(), *itr);
if (monitr != monitorSpecList.end())
++monCounter;
}
monitorwsSpecs += monCounter;
normalwsSpecs = m_total_specs - monitorwsSpecs;
g_log.debug() << "normalwsSpecs is " << normalwsSpecs << " monitorwsSpecs is "
<< monitorwsSpecs << std::endl;
Gigg, Martyn Anthony
committed
}
}//end if loop for m_interval
else
{ //if only List true
specid_t mSize = 0;
std::vector<specid_t>::const_iterator itr;
std::vector<specid_t>::const_iterator monitr;
Gigg, Martyn Anthony
committed
for (itr = m_spec_list.begin(); itr != m_spec_list.end(); ++itr)
{
monitr = find(monitorSpecList.begin(), monitorSpecList.end(), *itr);
if (monitr != monitorSpecList.end())
{
++mSize;
}
}
monitorwsSpecs = mSize;