Mantid
Loading...
Searching...
No Matches
LoadSQW.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
10#include "MantidAPI/Progress.h"
12#include "MantidAPI/Run.h"
13#include "MantidAPI/Sample.h"
22#include "MantidKernel/Matrix.h"
23#include "MantidKernel/Memory.h"
27#include <cfloat>
28
29using namespace Mantid::Kernel;
30using namespace Mantid::API;
32using namespace Mantid::DataObjects;
33
34namespace Mantid::MDAlgorithms {
35
36namespace {
46template <typename T> T interpretAs(std::vector<char> &Buf, size_t ind = 0) {
47 return *((reinterpret_cast<T *>(&Buf[ind])));
48}
49} // namespace
50
52
53
55 : m_fileName(""), m_fileStream(), m_prog(new Mantid::API::Progress(this, 0.05, 0.95, 100)), m_outputFile(""),
56 m_dataPositions(), m_boxSizes(), m_nDataPoints(0), m_mdImageSize(0), m_nDims(0), m_nBins() {}
57
65
66 // only .sqw can be considered
67 const std::string &extn = descriptor.extension();
68 if (extn != ".sqw")
69 return 0;
70
71 if (descriptor.isAscii()) {
72 return 10; // Low so that others may try
73 }
74 return 80; // probably it is sqw indeed
75}
76
79
81
84 std::vector<std::string> fileExtensions{".sqw"};
85 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, fileExtensions),
86 "File of type SQW format");
87 declareProperty(std::make_unique<API::WorkspaceProperty<API::IMDEventWorkspace>>("OutputWorkspace", "",
89 "Output IMDEventWorkspace reflecting SQW data read-in.");
90 declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("MetadataOnly", false),
91 "Load Metadata without events.");
92 std::vector<std::string> fileExtensions2{".nxs"};
94 std::make_unique<API::FileProperty>("OutputFilename", "", API::FileProperty::OptionalSave, fileExtensions2),
95 "If the input SQW file is too large to fit in memory, specify an output "
96 "NXS file.\n"
97 "The MDEventWorkspace will be create with this file as its back-end.");
98}
99
102
103 m_fileName = std::string(getProperty("Filename"));
104 // Parse Extract metadata. Including data locations.
106
107 // Create a new output workspace.
108 auto pWs = new MDEventWorkspace<MDEvent<4>, 4>;
110
111 // Add dimensions onto workspace.
112 std::vector<Mantid::Geometry::MDHistoDimensionBuilder> DimVector;
113
114 readDNDDimensions(DimVector, false);
115 readSQWDimensions(DimVector);
116 addDimsToWs(pWs, DimVector);
117 // Set some reasonable values for the box controller
118 BoxController_sptr bc = pWs->getBoxController();
119 for (size_t i = 0; i < 4; i++) {
120 bc->setSplitInto(i, m_nBins[i]);
121 }
122 bc->setMaxDepth(1);
123
124 // Initialize the workspace.
125 pWs->initialize();
126
127 // Add oriented lattice.
128 addLattice(pWs);
129
130 // Start with a MDGridBox.
131 pWs->splitBox();
132 //
133 readBoxSizes();
134
135 // Save the empty WS and turn it into a file-backed MDEventWorkspace (on
136 // option)
137 m_outputFile = getPropertyValue("OutputFilename");
138 if (!m_outputFile.empty()) {
139 // set file backed;
140 MemoryStats stat;
141 if ((m_nDataPoints * sizeof(MDEvent<4>) * 2 / 1024) < stat.availMem())
142 g_log.notice() << "You have enough memory available to load the " << m_nDataPoints
143 << " points into memory; this would be "
144 "faster than using a file back-end.\n";
145
146 auto saver = createChildAlgorithm("SaveMD", 0.01, 0.05, true);
147 saver->setProperty("InputWorkspace", ws);
148 saver->setPropertyValue("Filename", m_outputFile);
149 // should think about it.
150 saver->setProperty("UpdateFileBackEnd", false);
151 saver->setProperty("MakeFileBacked", false);
152 saver->executeAsChildAlg();
153 m_prog->resetNumSteps(100, 0.05, 0.75);
154
155 // set file backed boxes
156 auto Saver = std::shared_ptr<API::IBoxControllerIO>(new DataObjects::BoxControllerNeXusIO(bc.get()));
157 bc->setFileBacked(Saver, m_outputFile);
158 pWs->getBox()->setFileBacked();
159 bc->getFileIO()->setWriteBufferSize(1000000);
160
161 } else {
162 MemoryStats stat;
163 if ((size_t(double(m_nDataPoints) * 1.5) * sizeof(MDEvent<4>) / 1024) > stat.availMem())
164 g_log.warning() << "You may not have enough physical memory available to load the " << m_nDataPoints
165 << " points into memory. You can cancel and specify "
166 "OutputFilename to load to a file back-end.\n";
167 }
168
169 if (bc->isFileBacked()) {
170 std::cout << "File backed? " << bc->isFileBacked() << ". Cache " << bc->getFileIO()->getMemoryStr() << '\n';
171 } else {
172 bool ff(false);
173 std::cout << "File backed? " << ff << ". Cache 0\n";
174 }
175
176 // Persist the workspace.
177 API::IMDEventWorkspace_sptr i_out = getProperty("OutputWorkspace");
178 if (i_out)
179 throw std::runtime_error("Cannot currently handle overwriting of an existing workspace.");
180 setProperty("OutputWorkspace", ws);
181
182 // Read events into the workspace.
183 bool bMetadataOnly = getProperty("MetadataOnly");
184 if (!bMetadataOnly)
185 readEvents(pWs);
186
187 progress(m_outputFile.empty() ? 0.96 : 0.76, "Refreshing cache");
188 pWs->refreshCache();
189
190 if (!m_outputFile.empty()) {
191 g_log.notice() << "Starting SaveMD to update the file back-end.\n";
192 auto saver = createChildAlgorithm("SaveMD", 0.76, 1.00);
193 saver->setProperty("InputWorkspace", ws);
194 saver->setProperty("UpdateFileBackEnd", true);
195 saver->executeAsChildAlg();
196 }
197}
198
201 CPUTimer tim;
202
203 size_t maxNPix = ~size_t(0);
204 if (m_nDataPoints > maxNPix) {
205 throw std::runtime_error("Not possible to address all datapoints in "
206 "memory using this architecture ");
207 }
208
209 const size_t ncolumns = 9; // qx, qy, qz, en, s, err, irunid, idetid, ien
210 const size_t column_size = 4; // types stored as either float32 or unsigned integer (both 4byte).
211 const size_t column_size_2 = column_size * 2; // offset, gives qz
212 const size_t column_size_3 = column_size * 3; // offset, gives en
213 const size_t column_size_4 = column_size * 4; // offset, gives idet
214 // const size_t column_size_5 = column_size * 5; //offset, gives ien
215 const size_t column_size_6 = column_size * 6; // offset, gives irun
216 // TODO can this be read from the file? Was the file produced with SaveMD?
217 uint16_t goniometerIndex(0);
218 const size_t column_size_7 = column_size * 7; // offset, gives signal
219 const size_t column_size_8 = column_size * 8; // offset, gives error
220
221 const size_t pixel_width = ncolumns * column_size;
222 const size_t data_buffer_size = pixel_width * m_nDataPoints;
223 g_log.information() << m_nDataPoints << " data points in this SQW file.\n";
224
225 // Load from the input file is smallish blocks
226 size_t blockSize = pixel_width * 1000000;
227
228 // Report progress once per block
229 auto numBlocks = int((data_buffer_size + blockSize - 1) / blockSize);
230 m_prog->setNumSteps(numBlocks);
231 m_prog->setNotifyStep(0.1);
232
233 // For tracking when to split boxes
234 size_t eventsAdded = 0;
235 BoxController_sptr bc = ws->getBoxController();
236 DiskBuffer *dbuf(nullptr);
237 if (bc->isFileBacked())
238 dbuf = bc->getFileIO();
239
240 for (int blockNum = 0; blockNum < numBlocks; blockNum++) {
241 // Start point in the file
242 size_t inputFileOffset = size_t(blockNum) * blockSize;
243
244 // Limit the size of the block
245 size_t currentBlockSize = blockSize;
246 if ((inputFileOffset + currentBlockSize) > data_buffer_size)
247 currentBlockSize = data_buffer_size - inputFileOffset;
248
249 // Load the block from the file
250 std::vector<char> Buffer = std::vector<char>(currentBlockSize);
251 this->m_fileStream.seekg(this->m_dataPositions.pix_start + inputFileOffset, std::ios::beg);
252 this->m_fileStream.read(&Buffer[0], currentBlockSize);
253
254 // Go through each pixel in the input
255 auto currentNumPixels = int(currentBlockSize / pixel_width);
256 eventsAdded += size_t(currentNumPixels);
257
258 // Add the events in parallel
260 for (int i = 0; i < currentNumPixels; i++) {
261 auto current_pix = size_t(i * pixel_width);
262 coord_t centers[4] = {interpretAs<float>(Buffer, current_pix),
263 interpretAs<float>(Buffer, current_pix + column_size),
264 interpretAs<float>(Buffer, current_pix + column_size_2),
265 interpretAs<float>(Buffer, current_pix + column_size_3)};
266 const auto errorSQ = interpretAs<float>(Buffer, current_pix + column_size_8);
267 ws->addEvent(MDEvent<4>(
268 interpretAs<float>(Buffer, current_pix + column_size_7), // Signal
269 errorSQ, // Error sq
270 static_cast<uint16_t>(interpretAs<float>(Buffer, current_pix + column_size_6)), // run Index
271 goniometerIndex, static_cast<int32_t>(interpretAs<float>(Buffer, current_pix + column_size_4)), // Detector Id
272 centers));
273 }
274
275 // MemoryStats stat;
276 // size_t bytesAvail = stat.availMem() * 1024;
277 // // Estimate how many extra bytes will (temporarily) be used when
278 // splitting events
279 // size_t bytesNeededToSplit = eventsAdded * sizeof(MDEvent<4>) / 2;
280
281 // Split:
282 // 1. When < 1 GB of memory is free
283 // 2. When too little memory might be available for the splitting operation
284 // 3. When enough events have been added that it makes sense
285 // 4. At the last block being added
286
287 // if ((bytesAvail < 1000000000) || (bytesAvail < bytesNeededToSplit)
288 // ||
289 // bc->shouldSplitBoxes(eventsAdded*2, lastNumBoxes) || (blockNum
290 // == numBlocks-1) )
291
292 if (eventsAdded > 19000000) {
293 g_log.information() << "Splitting boxes after " << eventsAdded << " events added.\n";
294
295 // This splits up all the boxes according to split thresholds and sizes.
297 ThreadPool tp(ts);
298 ws->splitAllIfNeeded(ts);
299 tp.joinAll();
300
301 // Flush the cache - this will save things out to disk
302 if (dbuf)
303 dbuf->flushCache();
304 eventsAdded = 0;
305 }
306
307 // This will be correct optimized code after all.
308 //
309 // if (false)
310 // {
311 // std::vector<MDBoxBase<MDEvent<4>,4>*> boxes;
312 // ws->getBox()->getBoxes(boxes, 100, true);
313 // size_t modified = 0;
314 // size_t inmem = 0;
315 // size_t ondisk = 0;
316 // size_t events = 0;
317 // for (size_t i=0; i<boxes.size(); i++)
318 // {
319 // MDBox<MDEvent<4>,4>* box =
320 // dynamic_cast<MDBox<MDEvent<4>,4>*>(boxes[i]);
321 // if (box)
322 // {
323 // //box->save();
324 // if (box->dataAdded() || box->dataModified())
325 // modified++;
326 // if (box->getInMemory())
327 // inmem++;
328 // if (box->getOnDisk())
329 // ondisk++;
330 // events += box->getEventVectorSize();
331 // }
332 // }
333 // g_log.information() << modified << " of " << boxes.size() << "
334 // MDBoxes have data added or modified.\n";
335 // g_log.information() << inmem << " MDBoxes are in memory." <<
336 // '\n';
337 // //g_log.information() << ondisk << " MDBoxes are on disk." <<
338 // '\n';
339 // g_log.information() << double(events)/1e6 << " million events in
340 // memory.\n";
341 // }
342
343 // Report progress once per block.
344 m_prog->report();
345 }
347 ThreadPool tp(ts);
348 ws->splitAllIfNeeded(ts);
349 tp.joinAll();
350
351 // Flush the cache - this will save things out to disk
352 if (dbuf)
353 dbuf->flushCache();
354}
355
362 std::vector<char> buf(4 * (3 + 3)); // Where 4 = size_of(float) and 3 * 3 is size of b-matrix.
363 this->m_fileStream.seekg(this->m_dataPositions.geom_start, std::ios::beg);
364 this->m_fileStream.read(&buf[0], buf.size());
365
366 auto a = static_cast<double>(interpretAs<float>(buf, 0));
367 auto b = static_cast<double>(interpretAs<float>(buf, 4));
368 auto c = static_cast<double>(interpretAs<float>(buf, 8));
369 auto aa = static_cast<double>(interpretAs<float>(buf, 12));
370 auto bb = static_cast<double>(interpretAs<float>(buf, 16));
371 auto cc = static_cast<double>(interpretAs<float>(buf, 20));
372
374 // set up the goniometer. All mdEvents (pixels) in Horace sqw file are in lab
375 // frame,
376 // Q units so general goniometer should provide unit rotation matrix
377 info->mutableRun().mutableGoniometer().makeUniversalGoniometer();
378 info->mutableSample().setOrientedLattice(std::make_unique<OrientedLattice>(a, b, c, aa, bb, cc));
379 ws->addExperimentInfo(info);
380}
381
382void LoadSQW::buildMDDimsBase(std::vector<Mantid::Geometry::MDHistoDimensionBuilder> &DimVector) {
383 std::vector<std::string> dimID(4, "qx");
384 std::vector<std::string> dimUnit(4, "A^-1");
385 std::vector<std::string> dimFrameName(4, "HKL");
386 dimID[1] = "qy";
387 dimID[2] = "qz";
388 dimID[3] = "en";
389 dimUnit[3] = "meV";
390 dimFrameName[3] = "meV";
391 DimVector.resize(4);
392 for (size_t i = 0; i < 4; i++) {
393 DimVector[i].setId(dimID[i]);
394 DimVector[i].setUnits(dimUnit[i]);
395 DimVector[i].setName(dimID[i]);
396 DimVector[i].setFrameName(dimFrameName[i]);
397 }
398}
399
401void LoadSQW::readDNDDimensions(std::vector<Mantid::Geometry::MDHistoDimensionBuilder> &DimVectorOut,
402 bool arrangeByMDImage) {
403 // using Mantid::Geometry::MDHistoDimensionBuilder;
404 std::vector<Mantid::Geometry::MDHistoDimensionBuilder> DimVectorIn;
405 this->buildMDDimsBase(DimVectorIn);
406
407 std::vector<char> buf(4 * (3 + 3 + 4 + 16 + 4 + 2));
408 this->m_fileStream.seekg(this->m_dataPositions.geom_start, std::ios::beg);
409
410 this->m_fileStream.read(&buf[0], buf.size());
411 // skip allat and adlngldef
412 // interpret shifts
413 size_t i0 = 4 * (3 + 3);
414 // for(size_t i=0;i<this->m_nDims;i++){
415 // double val = (double)*((float*)(&buf[i0+i*4]));
416 // dscrptn.pDimDescription(i)->data_shift = val;
417 //}
418
419 // TODO: how to use it in our framework? -> it is B^-1 matrix possibly
420 // re-scaled
421 i0 += 4 * 4;
422 // [data.u_to_rlu, count, ok, mess] = fread_catch(fid,[4,4],'float32'); if
423 // ~all(ok); return; end;
424 size_t ic = 0;
425 for (size_t i = 0; i < 4; i++) {
426 for (size_t j = 0; j < 4; j++) {
427 ic++;
428 }
429 }
430 i0 += ic * 4;
431
432 // axis labels size
433 i0 += 4 * 4;
434 auto nRows = interpretAs<uint32_t>(buf, i0);
435 auto nCols = interpretAs<uint32_t>(buf, i0 + 4);
436
437 // read axis labelsg
438 buf.resize(nRows * nCols);
439 // [ulabel, count, ok, mess] = fread_catch(fid,[n(1),n(2)],'*char'); if
440 // ~all(ok); return; end;
441
442 this->m_fileStream.read(&buf[0], buf.size());
443
444 // data.ulabel=cellstr(ulabel)';
445 std::string name;
446 char symb;
447 name.resize(nCols);
448 for (unsigned int i = 0; i < nRows; i++) {
449 for (unsigned int j = 0; j < nCols; j++) {
450 symb = buf[i + j * nRows];
451 name[j] = symb;
452 }
453 // Trim string.
454 std::string sName(name);
455 boost::erase_all(sName, " ");
456
457 DimVectorIn[i].setName(sName);
458 }
459
460 // resize for iax and npax;
461 buf.resize(4 * 4 * 3);
462 this->m_fileStream.read(&buf[0], 4);
463
464 auto npax = interpretAs<uint32_t>(buf);
465 unsigned int niax = 4 - npax;
466
467 /*
468 [npax, count, ok, mess] = fread_catch(fid,1,'int32'); if ~all(ok); return;
469 end;
470 niax=4-npax;
471 if niax~=0
472 [data.iax, count, ok, mess] = fread_catch(fid,[1,niax],'int32'); if ~all(ok);
473 return; end;
474 [data.iint, count, ok, mess] = fread_catch(fid,[2,niax],'float32'); if
475 ~all(ok); return; end;
476 else
477 data.iax=zeros(1,0); % create empty index of integration array in standard
478 form
479 data.iint=zeros(2,0);
480 end
481 */
482 // axis counter
483 ic = 0;
484 std::vector<unsigned int> iax;
485 if (niax > 0) {
486 buf.resize(4 * (niax + 2 * niax));
487 iax.resize(niax);
488 this->m_fileStream.read(&buf[0], buf.size());
489
490 for (unsigned int i = 0; i < niax; i++) {
491 iax[i] = interpretAs<uint32_t>(buf, i * 4) - 1;
492 auto min = interpretAs<float>(buf, 4 * (niax + i * 2));
493 float max = interpretAs<float>(buf, 4 * (niax + i * 2 + 1)) * (1 + FLT_EPSILON);
494
495 DimVectorIn[ic].setNumBins(1);
496 DimVectorIn[ic].setMax(max);
497 DimVectorIn[ic].setMin(min);
498 ic++;
499 }
500 }
501
502 // processing projection axis;
503
504 /*
505 if npax~=0
506 [data.pax, count, ok, mess] = fread_catch(fid,[1,npax],'int32'); if ~all(ok);
507 return; end;
508 psize=zeros(1,npax); % will contain number of bins along each dimension of
509 plot axes
510 for i=1:npax
511 [np,count,ok,mess] = fread_catch(fid,1,'int32'); if ~all(ok); return; end;
512 [data.p{i},count,ok,mess] = fread_catch(fid,np,'float32'); if ~all(ok);
513 return; end;
514 psize(i)=np-1;
515 end
516 [data.dax, count, ok, mess] = fread_catch(fid,[1,npax],'int32'); if ~all(ok);
517 return; end;
518 if length(psize)==1
519 psize=[psize,1]; % make size of a column vector
520 end
521 else
522 data.pax=zeros(1,0); % create empty index of plot axes
523 data.p=cell(1,0);
524 data.dax=zeros(1,0); % create empty index of plot axes
525 psize=[1,1]; % to hold a scalar
526 end
527 */
528 std::vector<unsigned int> pax, dax;
529 if (npax > 0) {
530 //[data.pax, count, ok, mess] = fread_catch(fid,[1,npax],'int32'); if
531 //~all(ok); return; end;
532 this->m_fileStream.read(&buf[0], 4 * npax);
533 pax.resize(npax);
534 dax.resize(npax);
535 for (unsigned int i = 0; i < npax; i++) {
536 // projection axis indexes
537 pax[i] = interpretAs<uint32_t>(buf, 4 * i) - 1;
538
539 std::vector<char> axis_buffer(101 * 4);
540 this->m_fileStream.read(&axis_buffer[0], 4);
541 auto nAxisPoints = interpretAs<uint32_t>(axis_buffer);
542 if (axis_buffer.size() < nAxisPoints * 4)
543 axis_buffer.resize(nAxisPoints * 4);
544
545 this->m_fileStream.read(&axis_buffer[0], 4 * nAxisPoints);
546
547 auto min = interpretAs<float>(axis_buffer, 0);
548 float max = interpretAs<float>(axis_buffer, 4 * (nAxisPoints - 1)) * (1 + FLT_EPSILON);
549
550 DimVectorIn[ic].setNumBins(nAxisPoints - 1);
551 DimVectorIn[ic].setMax(max);
552 DimVectorIn[ic].setMin(min);
553 ic++;
554 }
555 //[data.dax, count, ok, mess] = fread_catch(fid,[1,npax],'int32'); if
556 //~all(ok); return; end;
557 this->m_fileStream.read(&buf[0], 4 * npax);
558 for (unsigned int i = 0; i < npax; i++)
559 dax[i] = interpretAs<uint32_t>(buf, 4 * i) - 1;
560 }
561 if (arrangeByMDImage) {
562
563 // Place dimensions to output vector in the correct dimensions order;
564 size_t dimIndex = 0;
565 DimVectorOut.resize(4);
566 for (size_t i = 0; i < npax; i++) {
567 DimVectorOut[dimIndex] = DimVectorIn[pax[dax[i]]];
568 dimIndex++;
569 }
570
571 for (size_t i = 0; i < niax; i++) {
572 DimVectorOut[dimIndex] = DimVectorIn[iax[i]];
573 dimIndex++;
574 }
575 } else // arrange according to sqw
576 {
577 DimVectorOut.assign(DimVectorIn.begin(), DimVectorIn.end());
578 }
579
580 // set up proper dimension bin numbers to use in further calculations
581 this->m_nBins.resize(4);
582 for (size_t i = 0; i < 4; i++) {
583 m_nBins[i] = DimVectorOut[i].getNumBins();
584 if (m_nBins[i] < 1)
585 m_nBins[i] = 1;
586 }
587}
590 std::vector<Mantid::Geometry::MDHistoDimensionBuilder> &DimVector) {
591 // Add dimensions to the workspace by invoking the dimension builders.
592 for (size_t i = 0; i < 4; i++) {
593 ws->addDimension(DimVector[i].create());
594 }
595}
596
598void LoadSQW::readSQWDimensions(std::vector<Mantid::Geometry::MDHistoDimensionBuilder> &DimVectorOut) {
600
601 this->buildMDDimsBase(DimVectorOut);
602
603 std::vector<char> buf(4 * (4 + 4));
604 this->m_fileStream.seekg(this->m_dataPositions.min_max_start, std::ios::beg);
605
606 this->m_fileStream.read(&buf[0], buf.size());
607
608 for (unsigned int i = 0; i < 4; i++) {
609 auto min = interpretAs<float>(buf, 4 * i * 2);
610 float max = interpretAs<float>(buf, 4 * (i * 2 + 1)) * (1 + FLT_EPSILON);
611 DimVectorOut[i].setNumBins(m_nBins[i]);
612 DimVectorOut[i].setMax(max);
613 DimVectorOut[i].setMin(min);
614 }
615}
616
617/*==================================================================================
618Region: Functions in the following region are candidates for refactoring. Copied
619from MD_FileHoraceReader
620==================================================================================*/
621
626void LoadSQW::parseMetadata(const std::string &fileName) {
627 if (m_fileStream.is_open())
628 m_fileStream.close();
629 m_fileStream.open(fileName.c_str(), std::ios::binary);
630
631 if (!m_fileStream.is_open())
632 throw(Kernel::Exception::FileError("Can not open input sqw file", fileName));
633
634 std::vector<char> data_buffer;
636 data_buffer.resize(3 * 4);
637
638 m_fileStream.read(&data_buffer[0], 2 * 4);
639 this->m_nDims = interpretAs<uint32_t>(data_buffer);
640
642
643 // go through all component headers and read them (or calculate their length)
644 std::streamoff next_position = m_dataPositions.component_headers_starts[0];
645 size_t nFiles = m_dataPositions.component_headers_starts.size();
646 for (size_t i = 0; i < nFiles; i++) {
647 m_dataPositions.component_headers_starts[i] = next_position;
648 next_position = m_dataPositions.parse_component_header(m_fileStream, next_position);
649 }
650 m_dataPositions.detectors_start = next_position;
651 // get detectors
653 // calculate all other data fields locations;
655}
656
658
660 m_fileStream.seekg(m_dataPositions.n_cell_pix_start, std::ios::beg);
661 m_fileStream.read(reinterpret_cast<char *>(&m_boxSizes[0]), m_dataPositions.mdImageSize * sizeof(uint64_t));
662}
663
664namespace LoadSQWHelper {
665
666// auxiliary functions
670void dataPositions::parse_sqw_main_header(std::ifstream &dataStream) { // we do not need this header at the moment ->
671 // just need to calculated its length;
672
673 std::vector<char> data_buffer(4 * 3);
674 dataStream.read(&data_buffer[0], 4);
675
676 unsigned int file_name_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
677 // skip main header file name
678 dataStream.seekg(file_name_length, std::ios_base::cur);
679
680 dataStream.read(&data_buffer[0], 4);
681 unsigned int file_path_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
682
683 // skip main header file path
684 dataStream.seekg(file_path_length, std::ios_base::cur);
685
686 dataStream.read(&data_buffer[0], 4);
687 unsigned int file_title = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
688
689 // skip ws title
690 dataStream.seekg(file_title, std::ios_base::cur);
691
692 // identify number of file headers, contributed into the dataset
693 dataStream.read(&data_buffer[0], 4);
694 unsigned int nFiles = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
695
697 this->component_headers_starts.assign(nFiles, 0);
698
699 std::streamoff last_location = dataStream.tellg();
700 if (last_location < 0)
701 throw("IO error for input file at start of component headers; Can not "
702 "seek to last location");
703 if (nFiles > 0) {
704 this->component_headers_starts[0] = last_location;
705 }
706}
714std::streamoff
716 std::streamoff start_location) { // we do not need this header at the
717 // moment -> just calculating its length;
718 // or may be we do soon?
719 std::vector<char> data_buffer(8);
720
721 std::streamoff shift = start_location - dataStream.tellg();
722 // move to specified location, which should be usually 0;
723 dataStream.seekg(shift, std::ios_base::cur);
724
725 dataStream.read(&data_buffer[0], 4);
726
727 unsigned int file_name_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
728 // skip component header file name
729 dataStream.seekg(file_name_length, std::ios_base::cur);
730
731 dataStream.read(&data_buffer[0], 4);
732 unsigned int file_path_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
733
734 // skip component header file path
735 dataStream.seekg(file_path_length, std::ios_base::cur);
736
737 // move to by specified number of bytes, see Matlab header above;
738 dataStream.seekg(4 * (7 + 3 * 4), std::ios_base::cur);
739
740 // read number of energy bins;
741 dataStream.read(&data_buffer[0], 4);
742 unsigned int nEn_bins = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
743 // skip energy values;
744 dataStream.seekg(4 * (nEn_bins), std::ios_base::cur);
745
746 // skip offsets and conversions;
747 dataStream.seekg(4 * (4 + 4 * 4 + 4), std::ios_base::cur);
748
749 // get labels matrix size;
750 dataStream.read(&data_buffer[0], 8);
751
752 unsigned int nRows = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
753 unsigned int nCols = *(reinterpret_cast<uint32_t *>(&data_buffer[4]));
754
755 // skip labels
756 dataStream.seekg(nRows * nCols, std::ios_base::cur);
757
758 std::streamoff end_location = dataStream.tellg();
759 return end_location;
760}
768std::streamoff dataPositions::parse_sqw_detpar(std::ifstream &dataStream,
769 std::streamoff start_location) { //
770 std::vector<char> data_buffer(8);
771
772 std::streamoff shift = start_location - dataStream.tellg();
773 // move to specified location, which should be usually 0;
774 dataStream.seekg(shift, std::ios_base::cur);
775
776 dataStream.read(&data_buffer[0], 4);
777 unsigned int file_name_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
778 // skip component header file name
779 dataStream.seekg(file_name_length, std::ios_base::cur);
780
781 dataStream.read(&data_buffer[0], 4);
782 unsigned int file_path_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
783 // skip component header file path
784 dataStream.seekg(file_path_length, std::ios_base::cur);
785
786 dataStream.read(&data_buffer[0], 4);
787 unsigned int num_detectors = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
788 // skip detector information
789 dataStream.seekg(num_detectors * 6 * 4, std::ios_base::cur);
790
791 std::streamoff end_location = dataStream.tellg();
792 return end_location;
793}
800void dataPositions::parse_data_locations(std::ifstream &dataStream, std::streamoff data_start,
801 std::vector<size_t> &nBins, uint64_t &nDataPoints) {
802 std::vector<char> data_buffer(12);
803
804 std::streamoff shift = data_start - dataStream.tellg();
805 // move to specified location, which should be usually 0;
806 dataStream.seekg(shift, std::ios_base::cur);
807
808 dataStream.read(&data_buffer[0], 4);
809
810 unsigned int file_name_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
811 // skip dummy file name
812 dataStream.seekg(file_name_length, std::ios_base::cur);
813
814 dataStream.read(&data_buffer[0], 4);
815 unsigned int file_path_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
816
817 // skip dummy file path
818 dataStream.seekg(file_path_length, std::ios_base::cur);
819
820 dataStream.read(&data_buffer[0], 4);
821 unsigned int data_title_length = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
822
823 // skip data title
824 dataStream.seekg(data_title_length, std::ios_base::cur);
825
826 this->geom_start = dataStream.tellg();
827
828 dataStream.seekg(4 * (3 + 3 + 4 + 16 + 4), std::ios_base::cur);
829
830 // get label information and skip labels;
831 dataStream.read(&data_buffer[0], 8);
832
833 unsigned int n_labels = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
834 unsigned int labels_length = *(reinterpret_cast<uint32_t *>(&data_buffer[4]));
835 dataStream.seekg(n_labels * labels_length, std::ios_base::cur);
836
837 this->npax_start = dataStream.tellg();
838
839 dataStream.read(&data_buffer[0], 4);
840
841 unsigned int npax = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
842 unsigned int niax = 4 - npax;
843 if (niax != 0) {
844 dataStream.seekg(3 * niax * 4, std::ios_base::cur);
845 }
846 if (npax != 0) {
847 nBins.resize(npax);
848
849 // skip projection axis
850 dataStream.seekg(npax * 4, std::ios_base::cur);
851
852 mdImageSize = 1;
853 for (unsigned int i = 0; i < npax; i++) {
854 dataStream.read(&data_buffer[0], 4);
855
856 unsigned int nAxisPoints = *(reinterpret_cast<uint32_t *>(&data_buffer[0]));
857 nBins[i] = nAxisPoints - 1;
858 mdImageSize *= nBins[i];
859 dataStream.seekg(nAxisPoints * 4, std::ios_base::cur);
860 }
861 // skip display indexes;
862 dataStream.seekg(npax * 4, std::ios_base::cur);
863 }
864 // signal start:
865 this->s_start = dataStream.tellg();
866 // and skip to errors
867 dataStream.seekg(mdImageSize * 4, std::ios_base::cur);
868
869 // error start:
870 this->err_start = dataStream.tellg();
871 dataStream.seekg(mdImageSize * 4, std::ios_base::cur);
872
873 // dnd data file. we do not support this?
874 if (dataStream.eof()) {
875 nDataPoints = 0;
876 return;
877 // throw(std::invalid_argument("DND Horace datasets are not supported by
878 // Mantid"));
879 }
880
881 this->n_cell_pix_start = dataStream.tellg();
882 // skip to the end of pixels;
883 dataStream.seekg(mdImageSize * 8, std::ios_base::cur);
884
885 if (dataStream.eof()) {
886 nDataPoints = 0;
887 return;
888 // throw(std::invalid_argument("DND b+ Horace datasets are not supported by
889 // Mantid"));
890 }
891 this->min_max_start = dataStream.tellg();
892 // skip min-max start
893 //[data.urange,count,ok,mess] = fread_catch(fid,[2,4],'float32'); if ~all(ok);
894 // return; end;
895 dataStream.seekg(8 * 4, std::ios_base::cur);
896
897 if (dataStream.eof()) {
898 nDataPoints = 0;
899 return;
900 // throw(std::invalid_argument("SQW a- Horace datasets are not supported by
901 // Mantid"));
902 }
903 // skip redundant field and read nPix (number of data points)
904 dataStream.read(&data_buffer[0], 12);
905
906 nDataPoints = static_cast<size_t>(*(reinterpret_cast<uint64_t *>(&data_buffer[4])));
907 this->pix_start = dataStream.tellg();
908}
909
910/*==================================================================================
911EndRegion:
912==================================================================================*/
913} // namespace LoadSQWHelper
914} // namespace Mantid::MDAlgorithms
static std::unique_ptr< QThreadPool > tp
#define PARALLEL_FOR_NO_WSP_CHECK()
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
This class is shared by a few Workspace types and holds information related to a particular experimen...
@ OptionalSave
to specify a file to write to but an empty string is
Definition: FileProperty.h:50
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
The class responsible for saving events into nexus file using generic box controller interface Expect...
Templated class for the multi-dimensional event workspace.
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Definition: MDEvent.h:36
Class to implement UB matrix.
CPUTimer : Timer that uses the CPU time, rather than wall-clock time to measure execution time.
Definition: CPUTimer.h:24
Buffer objects that need to be written out to disk so as to optimize writing operations.
Definition: DiskBuffer.h:42
void flushCache()
Flush out all the data in the memory; and writes out everything in the to-write cache.
Definition: DiskBuffer.cpp:192
Records the filename and the description of failure.
Definition: Exception.h:98
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.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
This class is responsible for memory statistics.
Definition: Memory.h:28
std::size_t availMem() const
Returns the available memory of the system in kiB.
Definition: Memory.cpp:398
void resetNumSteps(int64_t nsteps, double start, double end)
Change the number of steps between start/end.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
void setNumSteps(int64_t nsteps)
Change the number of steps between start/end.
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,...
Definition: ThreadPool.h:36
A First-In-First-Out Thread Scheduler.
The ThreadScheduler object defines how tasks are allocated to threads and in what order.
void parseMetadata(const std::string &fileName)
Parse metadata from file.
Definition: LoadSQW.cpp:626
std::ifstream m_fileStream
File stream containing binary file data.
Definition: LoadSQW.h:136
std::string m_outputFile
OutputFilename param.
Definition: LoadSQW.h:142
~LoadSQW() override
Destructor.
Definition: LoadSQW.cpp:78
void readDNDDimensions(std::vector< Mantid::Geometry::MDHistoDimensionBuilder > &DimVectorOut, bool arrangeByMDImage=true)
Read DND dimensions.
Definition: LoadSQW.cpp:401
LoadSQWHelper::dataPositions m_dataPositions
Instance of helper type, which describes the positions of the data within binary Horace file.
Definition: LoadSQW.h:146
void addDimsToWs(Mantid::DataObjects::MDEventWorkspace< DataObjects::MDEvent< 4 >, 4 > *ws, std::vector< Mantid::Geometry::MDHistoDimensionBuilder > &DimVector)
add range of dimensions to the workspace;
Definition: LoadSQW.cpp:589
void readSQWDimensions(std::vector< Mantid::Geometry::MDHistoDimensionBuilder > &DimVectorOut)
Read SQW dimensions.
Definition: LoadSQW.cpp:598
std::vector< uint64_t > m_boxSizes
Definition: LoadSQW.h:148
std::string m_fileName
the name of the file to work with
Definition: LoadSQW.h:134
void init() override
Provide wiki documentation.
Definition: LoadSQW.cpp:83
std::vector< size_t > m_nBins
number of bins in every non-integrated dimension
Definition: LoadSQW.h:154
const std::string name() const override
Algorithm's name for identification.
Definition: LoadSQW.h:89
virtual void readEvents(Mantid::DataObjects::MDEventWorkspace< DataObjects::MDEvent< 4 >, 4 > *ws)
Read events onto the workspace.
Definition: LoadSQW.cpp:200
Mantid::API::Progress * m_prog
Progress bar.
Definition: LoadSQW.h:139
void readBoxSizes()
read real box sizes for all Horace Boxes;
Definition: LoadSQW.cpp:657
void exec() override
Execute the algorithm.
Definition: LoadSQW.cpp:101
int confidence(Kernel::FileDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
Definition: LoadSQW.cpp:64
virtual void addLattice(Mantid::DataObjects::MDEventWorkspace< DataObjects::MDEvent< 4 >, 4 > *ws)
Extract lattice information.
Definition: LoadSQW.cpp:361
void buildMDDimsBase(std::vector< Mantid::Geometry::MDHistoDimensionBuilder > &DimVector)
build an initial range of 4 dimensions
Definition: LoadSQW.cpp:382
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< BoxController > BoxController_sptr
Shared ptr to BoxController.
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
Helper class which provides the Collimation Length for SANS instruments.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
@ Output
An output workspace.
Definition: Property.h:54
std::vector< std::streamoff > component_headers_starts
Definition: LoadSQW.h:41
std::streamoff parse_sqw_detpar(std::ifstream &dataStream, std::streamoff start_location)
Block 3: Detpar: parse positions of the contributed detectors.
Definition: LoadSQW.cpp:768
std::streamoff parse_component_header(std::ifstream &dataStream, std::streamoff start_location)
Block 2: Header: Parse header of single SPE file.
Definition: LoadSQW.cpp:715
void parse_sqw_main_header(std::ifstream &dataStream)
Block 1: Main_header: Parse SQW main data header.
Definition: LoadSQW.cpp:670
void parse_data_locations(std::ifstream &dataStream, std::streamoff data_start, std::vector< size_t > &nBins, uint64_t &nDataPoints)
Block 4: Data: parse positions of the data fields.
Definition: LoadSQW.cpp:800