30using API::ExperimentInfo;
31using Geometry::Goniometer;
32using Geometry::MDHistoDimensionBuilder;
33using Geometry::OrientedLattice;
34using Kernel::BinaryStreamReader;
43constexpr int64_t NPIX_CHUNK = 150000;
46constexpr int64_t NCHUNKS_SPLIT = 125;
48constexpr int32_t FIELDS_PER_PIXEL = 9;
50constexpr double INV_TWO_PI = 0.5 / M_PI;
63const std::string
LoadSQW2::category()
const {
return "DataHandling\\SQW;MDAlgorithms\\DataHandling"; }
68 return "Load an N-dimensional workspace from a .sqw file produced by "
80 const std::string &extn = descriptor.
extension();
97 using StringInitializerList = std::initializer_list<std::string>;
101 "File of type SQW format");
105 "If specified, the output workspace will be a file-backed "
107 std::vector<std::string> allowed = {
"Q_sample",
"HKL"};
108 declareProperty(
"Q3DFrames", allowed[0], std::make_shared<StringListValidator>(allowed),
109 "The required frame for the output workspace");
114 "Output IMDEventWorkspace reflecting SQW data");
148 std::string appName, filename, filepath, title;
149 double appVersion(0.0);
150 int32_t sqwType(-1), numDims(-1), nspe(-1);
151 *
m_reader >> appName >> appVersion >> sqwType >> numDims >> filename >> filepath >> title >> nspe;
152 m_nspe =
static_cast<uint16_t
>(nspe);
153 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
154 std::ostringstream os;
155 os <<
"Main header:\n"
156 <<
" app_name: " << appName <<
"\n"
157 <<
" app_version: " << appVersion <<
"\n"
158 <<
" sqw_type: " << sqwType <<
"\n"
159 <<
" ndims: " << numDims <<
"\n"
160 <<
" filename: " << filename <<
"\n"
161 <<
" filepath: " << filepath <<
"\n"
162 <<
" title: " << title <<
"\n"
163 <<
" nfiles: " <<
m_nspe <<
"\n";
175 throw std::runtime_error(
"Unsupported SQW type: " +
std::to_string(sqwType) +
176 "\nOnly files containing the full pixel "
177 "information are currently supported");
191 for (uint16_t i = 0; i <
m_nspe; ++i) {
195 auto expt0 =
m_outputWS->getExperimentInfo(0);
206 auto experiment = std::make_shared<ExperimentInfo>();
207 auto &sample = experiment->mutableSample();
208 auto &run = experiment->mutableRun();
217 run.addProperty(
"Ei",
static_cast<double>(efix),
true);
220 std::vector<float> floats;
222 auto lattice = std::make_unique<OrientedLattice>(floats[0], floats[1], floats[2], floats[3], floats[4], floats[5]);
223 V3D uVec(floats[6], floats[7], floats[8]), vVec(floats[9], floats[10], floats[11]);
224 lattice->setUFromVectors(uVec, vVec);
225 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
226 std::stringstream os;
228 <<
" alatt: " << lattice->a1() <<
" " << lattice->a2() <<
" " << lattice->a3() <<
"\n"
229 <<
" angdeg: " << lattice->alpha() <<
" " << lattice->beta() <<
" " << lattice->gamma() <<
"\n"
230 <<
" cu: " << floats[6] <<
" " << floats[7] <<
" " << floats[8] <<
"\n"
231 <<
" cv: " << floats[9] <<
" " << floats[10] <<
" " << floats[11] <<
"\n"
232 <<
"B matrix (calculated): " << lattice->getB() <<
"\n"
233 <<
"Inverse B matrix (calculated): " << lattice->getBinv() <<
"\n";
236 sample.setOrientedLattice(std::move(lattice));
239 float psi(0.0f), omega(0.0f), dpsi(0.0f), gl(0.0f), gs(0.0f);
240 *
m_reader >> psi >> omega >> dpsi >> gl >> gs;
243 goniometer.
pushAxis(
"psi", uvCross[0], uvCross[1], uvCross[2], psi);
244 goniometer.
pushAxis(
"omega", uvCross[0], uvCross[1], uvCross[2], omega);
245 goniometer.
pushAxis(
"gl", 1.0, 0.0, 0.0, gl);
246 goniometer.
pushAxis(
"gs", 0.0, 0.0, 1.0, gs);
247 goniometer.
pushAxis(
"dpsi", 0.0, 1.0, 0.0, dpsi);
248 run.setGoniometer(goniometer,
false);
249 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
250 std::stringstream os;
251 os <<
"Goniometer angles:\n"
252 <<
" psi: " << psi <<
"\n"
253 <<
" omega: " << omega <<
"\n"
254 <<
" gl: " << gl <<
"\n"
255 <<
" gs: " << gs <<
"\n"
256 <<
" dpsi: " << dpsi <<
"\n"
257 <<
" goniometer matrix: " << goniometer.
getR() <<
"\n";
263 std::vector<float> enBins(nbounds);
265 run.storeHistogramBinBoundaries(std::vector<double>(enBins.begin(), enBins.end()));
269 m_file->seekg(96, std::ios_base::cur);
270 std::vector<int32_t> ulabel_shape(2);
273 m_file->seekg(ulabel_shape[0] * ulabel_shape[1], std::ios_base::cur);
291 std::string filename, filepath;
295 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
296 std::stringstream os;
297 os <<
"Skipping " << ndet <<
" detector parameters from '" << filename <<
"'\n";
301 m_file->seekg(6 * 4 * ndet, std::ios_base::cur);
319 *
m_reader >> dropped >> dropped >> dropped;
321 m_file->seekg(120, std::ios_base::cur);
324 std::vector<int32_t> ulabelShape(2);
326 m_file->seekg(ulabelShape[0] * ulabelShape[1], std::ios_base::cur);
337 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
338 std::stringstream os;
340 for (
const auto &val : nbins) {
347 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
348 std::stringstream os;
349 os <<
"data extents (in output frame): ";
350 for (
size_t i = 0; i < 4; ++i) {
351 os <<
"(" << dimLimits[2 * i] <<
"," << dimLimits[2 * i + 1] <<
") ";
360 const auto &bmat0 =
m_outputWS->getExperimentInfo(0)->sample().getOrientedLattice().getB();
361 for (
size_t i = 0; i < 4; ++i) {
364 float umin(dimLimits[2 * i]), umax(dimLimits[2 * i + 1]);
381 int32_t nProjAxes(0);
383 int32_t nIntAxes(4 - nProjAxes);
386 m_file->seekg(nIntAxes *
sizeof(int32_t) + 2 * nIntAxes *
sizeof(
float), std::ios_base::cur);
388 std::vector<int32_t> nbins(4, 1);
391 std::vector<int32_t> projAxIdx;
392 int32_t signalLength(1);
393 m_reader->read(projAxIdx, nProjAxes);
394 for (int32_t i = 0; i < nProjAxes; ++i) {
397 m_file->seekg(nbounds *
sizeof(
float), std::ios_base::cur);
398 nbins[projAxIdx[i] - 1] = nbounds - 1;
399 signalLength *= nbounds - 1;
402 m_file->seekg(nProjAxes *
sizeof(int32_t), std::ios_base::cur);
404 m_file->seekg(2 * signalLength *
sizeof(
float) + signalLength *
sizeof(int64_t), std::ios_base::cur);
421 m_file->seekg(8 *
sizeof(
float), std::ios_base::cur);
423 auto filePosAfterURange =
m_file->tellg();
425 m_file->seekg(
sizeof(int32_t), std::ios_base::cur);
432 constexpr int64_t bufferSize(FIELDS_PER_PIXEL * NPIX_CHUNK);
433 std::vector<float> pixBuffer(bufferSize);
434 int64_t pixelsLeftToRead(npixtot);
435 std::vector<float> dimLimits(8);
436 dimLimits[0] = dimLimits[2] = dimLimits[4] = dimLimits[6] = FLT_MAX;
437 dimLimits[1] = dimLimits[3] = dimLimits[5] = dimLimits[7] = -FLT_MAX;
438 while (pixelsLeftToRead > 0) {
439 int64_t chunkSize(pixelsLeftToRead);
440 if (chunkSize > NPIX_CHUNK) {
441 chunkSize = NPIX_CHUNK;
443 m_reader->read(pixBuffer, FIELDS_PER_PIXEL * chunkSize);
444 for (int64_t i = 0; i < chunkSize; ++i) {
445 float *pixel = pixBuffer.data() + i * 9;
447 for (
size_t j = 0; j < 4; ++j) {
449 if (uj < dimLimits[2 * j])
450 dimLimits[2 * j] = uj;
451 else if (uj > dimLimits[2 * j + 1])
452 dimLimits[2 * j + 1] = uj;
454 status.
report(
"Calculating data extents");
456 pixelsLeftToRead -= chunkSize;
458 m_file->seekg(filePosAfterURange);
479 throw std::logic_error(
"LoadSQW2::createQDimension - Expected a dimension "
480 "index between 0 & 2. Found: " +
483 static std::array<const char *, 3> indexToDim{
"x",
"y",
"z"};
485 builder.
setId(std::string(
"q") + indexToDim[
index]);
491 std::string name, unit, frameName;
492 if (m_outputFrame ==
"Q_sample") {
493 name = m_outputFrame +
"_" + indexToDim[
index];
495 frameName =
"QSample";
496 }
else if (m_outputFrame ==
"HKL") {
497 static std::array<const char *, 3> indexToHKL{
"[H,0,0]",
"[0,K,0]",
"[0,0,L]"};
498 name = indexToHKL[
index];
501 const V3D x = bmat * dimDir;
502 double length = 2. * M_PI *
x.norm();
506 throw std::logic_error(
"LoadSQW2::createQDimension - Unknown output frame: " + m_outputFrame);
544 auto boxController =
m_outputWS->getBoxController();
545 for (
size_t i = 0; i < 4; i++) {
546 boxController->setSplitInto(i,
m_outputWS->getDimension(i)->getNBins());
548 boxController->setMaxDepth(1);
555 std::string fileback =
getProperty(
"OutputFilename");
556 if (!fileback.empty()) {
569 savemd->setProperty(
"InputWorkspace",
m_outputWS);
570 savemd->setPropertyValue(
"Filename", filebackPath);
571 savemd->setProperty(
"UpdateFileBackEnd",
false);
572 savemd->setProperty(
"MakeFileBacked",
false);
573 savemd->executeAsChildAlg();
576 auto boxControllerMem =
m_outputWS->getBoxController();
577 auto boxControllerIO = std::make_shared<BoxControllerNeXusIO>(boxControllerMem.get());
578 boxControllerMem->setFileBacked(boxControllerIO, filebackPath);
580 boxControllerMem->getFileIO()->setWriteBufferSize(1000000);
591 m_file->seekg(
sizeof(int32_t), std::ios_base::cur);
594 g_log.
debug() <<
" npixtot: " << npixtot <<
"\n";
603 constexpr int64_t bufferSize(FIELDS_PER_PIXEL * NPIX_CHUNK);
604 std::vector<float> pixBuffer(bufferSize);
605 int64_t pixelsLeftToRead(npixtot), chunksRead(0);
606 size_t pixelsAdded(0);
607 while (pixelsLeftToRead > 0) {
608 int64_t chunkSize(pixelsLeftToRead);
609 if (chunkSize > NPIX_CHUNK) {
610 chunkSize = NPIX_CHUNK;
612 m_reader->read(pixBuffer, FIELDS_PER_PIXEL * chunkSize);
613 for (int64_t i = 0; i < chunkSize; ++i) {
615 status.
report(
"Reading pixel data to workspace");
617 pixelsLeftToRead -= chunkSize;
619 if ((chunksRead % NCHUNKS_SPLIT) == 0) {
623 assert(pixelsLeftToRead == 0);
624 if (pixelsAdded == 0) {
625 throw std::runtime_error(
"No pixels could be added from the source file. "
626 "Please check the irun fields of all pixels are valid.");
627 }
else if (pixelsAdded !=
static_cast<size_t>(npixtot)) {
628 g_log.
warning(
"Some pixels within the source file had an invalid irun "
629 "field. They have been ignored.");
660 size_t reqdMemory = (npixtot *
sizeof(
MDEvent<4>) + NPIX_CHUNK * FIELDS_PER_PIXEL) / 1024;
662 g_log.
warning() <<
"It looks as if there is insufficient memory to load the "
663 <<
"entire file. It is recommended to cancel the algorithm and "
665 "the OutputFilename option to create a file-backed workspace.\n";
679 uint16_t goniometerIndex(0);
683 auto irun =
static_cast<uint16_t
>(pixel[4]);
684 if (irun < 1 || irun >
m_nspe) {
687 coord_t centers[4] = {pixel[0], pixel[1], pixel[2], pixel[3]};
689 auto error = pixel[8];
691 goniometerIndex,
static_cast<detid_t>(pixel[5]), centers));
709 centers[0] =
static_cast<float>(qout[0]);
710 centers[1] =
static_cast<float>(qout[1]);
711 centers[2] =
static_cast<float>(qout[2]);
723 savemd->setProperty(
"InputWorkspace",
m_outputWS);
724 savemd->setProperty(
"UpdateFileBackEnd",
true);
725 savemd->executeAsChildAlg();
static std::unique_ptr< QThreadPool > tp
std::map< DeltaEMode::Type, std::string > index
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
#define GNU_DIAG_OFF(x)
This is a collection of macros for turning compiler warnings off in a controlled manner.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
@ OptionalSave
to specify a file to write to but an empty string is
@ Load
allowed here which will be passed to the algorithm
Helper class for reporting progress from algorithms.
A property class for workspaces.
The class responsible for saving events into nexus file using generic box controller interface Expect...
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Class to represent a particular goniometer setting, which is described by the rotation matrix.
void pushAxis(const std::string &name, double axisx, double axisy, double axisz, double angle=0., int sense=CCW, int angUnit=angDegrees)
Add an additional axis to the goniometer, closer to the sample.
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
MDHistoDimensionBuilder :
void setId(std::string id)
void setFrameName(std::string frameName)
Setter for the frame name.
IMDDimension_sptr create()
void setNumBins(size_t nbins)
void setName(const std::string &name)
static void resizeToFitMDBox(CoordT &min, CoordT &max)
Push the min/max values out by a defined amount.
void setUnits(const Kernel::UnitLabel &units)
Class to implement UB matrix.
const Kernel::DblMatrix & getBinv() const
Get the inverse of the B-matrix.
Defines a wrapper around an open file.
static bool isAscii(const std::string &filename, const size_t nbytes=256)
Returns true if the file is considered ascii.
const std::string & extension() const
Access the file extension.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
void debug(const std::string &msg)
Logs at debug level.
void warning(const std::string &msg)
Logs at warning level.
bool is(int level) const
Returns true if at least the given log level is set.
This class is responsible for memory statistics.
std::size_t availMem() const
Returns the available memory of the system in kiB.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
void setNotifyStep(double notifyStepPct)
Override the frequency at which notifications are sent out.
The concrete, templated class for properties.
A Thread Pool implementation that keeps a certain number of threads running (normally,...
A First-In-First-Out Thread Scheduler.
A simple class that provides a wall-clock (not processor time) timer.
float elapsed(bool reset=true)
Returns the wall-clock time elapsed in seconds since the Timer object's creation, or the last call to...
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
void setupFileBackend(const std::string &filebackPath)
Setup the filebackend for the output workspace.
void exec() override
Execute the algorithm.
void skipDetectorSection()
Skip the data in the detector section.
int32_t readMainHeader()
Reads the initial header section.
std::unique_ptr< Kernel::BinaryStreamReader > m_reader
std::shared_ptr< API::ExperimentInfo > readSingleSPEHeader()
Read single SPE header from the file.
void initFileReader()
Opens the file given to the algorithm and initializes the reader.
void readSQWDimensions()
Read and create the SQW dimensions on the output.
size_t addEventFromBuffer(const float *pixel)
Assume the given pointer points to the start of a full pixel and create an MDEvent based on it iff it...
std::unique_ptr< std::ifstream > m_file
void setupBoxController()
Setup the box controller based on the bin structure.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Kernel::DblMatrix m_uToRLU
Geometry::IMDDimension_sptr createEnDimension(float umin, float umax, size_t nbins)
Create an energy dimension.
void finalize()
Assumed to be the last step in the algorithm.
int version() const override
Algorithm's version for identification.
void init() override
Initialize the algorithm's properties.
std::vector< float > calculateDimLimitsFromData()
Find the dimension limits for each dimension in the target frame.
void skipDataSectionMetadata()
Skip metadata in data section.
Geometry::IMDDimension_sptr createQDimension(size_t index, float dimMin, float dimMax, size_t nbins, const Kernel::DblMatrix &bmat)
Create the Q MDHistoDimension for the output frame and given information from the file.
void warnIfMemoryInsufficient(int64_t npixtot)
If the output is not file backed and the machine appears to have insufficient memory to read the data...
void throwIfUnsupportedFileType(int32_t sqwType)
Throw std::runtime_error if the sqw type of the file is unsupported.
void cacheFrameTransforms(const Geometry::OrientedLattice &lattice)
Cache the transforms between the Q_sample & HKL frames from the given lattice.
void readAllSPEHeadersToWorkspace()
Read all of the SPE headers and fill in the experiment details on the output workspace.
void createOutputWorkspace()
Create the output workspace object.
void readPixelDataIntoWorkspace()
Read the pixel data into the workspace.
int confidence(Kernel::FileDescriptor &descriptor) const override
Return the confidence with this algorithm can load the file.
std::shared_ptr< SQWWorkspace > m_outputWS
void toOutputFrame(coord_t *centers)
Transform the given coordinates to the requested output frame if necessary.
void splitAllBoxes()
Split boxes in the output workspace if required.
const std::string category() const override
Algorithm's category for identification.
void cacheInputs()
Cache any user input to avoid repeated lookups.
std::vector< int32_t > readProjection()
Read the required parts of the projection information from the data section The file pointer is assum...
std::string m_outputFrame
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IMDDimension > IMDDimension_sptr
Shared Pointer for IMDDimension. Frequently used type in framework.
Mantid::Kernel::Matrix< double > DblMatrix
std::string DLLExport sprintfd(const double data, const double eps)
creates string representation of the number with accuracy, cpecified by eps
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
int32_t detid_t
Typedef for a detector ID.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Output
An output workspace.