Mantid
Loading...
Searching...
No Matches
LoadDNSEvent.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2022 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 +
7
12#include "MantidAPI/Run.h"
23#include "MantidKernel/System.h"
24#include <stdexcept>
25#include <vector>
26
27#include <chrono>
28#include <iostream>
29
30using separator_t = std::array<uint8_t, 8>;
31static constexpr separator_t header_sep{0x00, 0x00, 0x55, 0x55, 0xAA, 0xAA, 0xFF, 0xFF};
32static constexpr separator_t block_sep = {0x00, 0x00, 0xFF, 0xFF, 0x55, 0x55, 0xAA, 0xAA}; // 0xAAAA5555FFFF0000; //
33static constexpr separator_t closing_sig = {0xFF, 0xFF, 0xAA, 0xAA, 0x55, 0x55, 0x00, 0x00}; // 0x00005555AAAAFFFF; //
34
35using namespace Mantid::DataObjects;
36using namespace Mantid::Kernel;
37using namespace Mantid::API;
38
39namespace {
40
41template <typename V1, typename V2> bool startsWith(const V1 &sequence, const V2 &subSequence) {
42 return std::equal(std::begin(subSequence), std::end(subSequence), std::begin(sequence));
43}
44
45} // namespace
46
47namespace Mantid::DataHandling {
48
50
51
57int LoadDNSEvent::confidence(Kernel::FileDescriptor &descriptor) const {
58 const std::string &extn = descriptor.extension();
59 if (extn != ".mdat")
60 return 0;
61 auto &file = descriptor.data();
62
63 std::string fileline;
64 std::getline(file, fileline);
65 if (fileline.find("mesytec psd listmode data") != std::string::npos) {
66 return 80;
67 } else
68 return 0;
69}
70
71const unsigned int MAX_BUFFER_BYTES_SIZE = 1500; // maximum buffer size in data file
72const unsigned int PIXEL_PER_TUBE = 1024; // maximum buffer size in data file
73
76
77 const std::vector<std::string> exts{".mdat"};
78 declareProperty(std::make_unique<FileProperty>("InputFile", "", FileProperty::Load, exts),
79 "A DNS mesydaq listmode event datafile.");
80
81 declareProperty<uint32_t>("ChopperChannel", 2u, std::make_shared<BoundedValidator<uint32_t>>(1, 4),
82 "The Chopper Channel (1 to 4)", Kernel::Direction::Input);
83
84 declareProperty<uint32_t>("NumberOfTubes", 128, std::make_shared<BoundedValidator<uint32_t>>(1, 128),
85 "The number of tubes, each tube has 1024 pixels (1 to 128)", Kernel::Direction::Input);
86
87 declareProperty<bool>("DiscardPreChopperEvents", true, std::make_shared<BoundedValidator<bool>>(0, 1),
88 "Discards events before first chopper trigger (turn off for elastic)",
90
91 declareProperty<bool>("SetBinBoundary", true, std::make_shared<BoundedValidator<bool>>(0, 1),
92 "Sets all bin boundaries to include all events (can be turned off to save time).",
94
96 std::make_unique<WorkspaceProperty<DataObjects::EventWorkspace>>("OutputWorkspace", "", Direction::Output),
97 "The name of the output workspace.");
98}
99
102 // loadProperties:
103 const std::string fileName = getPropertyValue("InputFile");
104 m_chopperChannel = static_cast<uint32_t>(getProperty("ChopperChannel"));
105 m_discardPreChopperEvents = getProperty("DiscardPreChopperEvents");
106 m_setBinBoundary = getProperty("SetBinBoundary");
107 m_detectorPixelCount = static_cast<uint32_t>(getProperty("NumberOfTubes")) * PIXEL_PER_TUBE;
108
109 // create workspace
110 EventWorkspace_sptr outputWS = std::dynamic_pointer_cast<EventWorkspace>(
111 WorkspaceFactory::Instance().create("EventWorkspace", m_detectorPixelCount, 1, 1));
112 outputWS->switchEventType(Mantid::API::EventType::TOF);
113 outputWS->getAxis(0)->setUnit("TOF");
114 outputWS->setYUnit("Counts");
115
116 // g_log.notice() << "ChopperChannel: " << m_chopperChannel << std::endl;
117
118 FileByteStream file(static_cast<std::string>(fileName));
119
120 auto finalEventAccumulator = parse_File(file, fileName);
121 populate_EventWorkspace(outputWS, finalEventAccumulator);
122 if (m_setBinBoundary && (outputWS->getNumberEvents() > 0)) {
123 outputWS->setAllX({0, outputWS->getEventXMax()});
124 }
125 setProperty("OutputWorkspace", outputWS);
126}
127
129 static const unsigned EVENTS_PER_PROGRESS = 100;
130 // The number of steps depends on the type of input file
131 Progress progress(this, 0.0, 1.0, finalEventAccumulator.neutronEvents.size() / EVENTS_PER_PROGRESS);
132
133 // Sort reversed (latest event first, most early event last):
134 std::sort(finalEventAccumulator.triggerEvents.begin(), finalEventAccumulator.triggerEvents.end(),
135 [](auto l, auto r) { return l.timestamp > r.timestamp; });
136
137 std::atomic<unsigned int> oversizedChanelIndexCounterA(0);
138 std::atomic<unsigned int> oversizedPosCounterA(0);
139
141 for (int eventIndex = 0; eventIndex < static_cast<int>(finalEventAccumulator.neutronEvents.size()); ++eventIndex) {
142 // uint64_t chopperTimestamp = 0;
143 unsigned int oversizedChanelIndexCounter = 0;
144 unsigned int oversizedPosCounter = 0;
145 const auto wsIndex = eventIndex;
146 auto &eventList = finalEventAccumulator.neutronEvents[eventIndex];
147 if (eventList.size() != 0) {
148 std::sort(eventList.begin(), eventList.end(), [](auto l, auto r) { return l.timestamp < r.timestamp; });
149 }
150
151 auto chopperIt = finalEventAccumulator.triggerEvents.cbegin();
152
153 auto &spectrum = eventWS->getSpectrum(wsIndex);
155
156 uint64_t numProcessed = 0;
157 for (const auto &event : eventList) {
158 numProcessed++;
159 if (numProcessed % EVENTS_PER_PROGRESS == 0) {
160 progress.report();
161 if (this->getCancel()) {
162 throw CancelException();
163 }
164 }
165
166 chopperIt =
167 std::lower_bound(finalEventAccumulator.triggerEvents.cbegin(), finalEventAccumulator.triggerEvents.cend(),
168 event.timestamp, [](auto l, auto r) { return l.timestamp > r; });
169 const uint64_t chopperTimestamp = chopperIt != finalEventAccumulator.triggerEvents.cend()
170 ? chopperIt->timestamp
171 : 0; // before first chopper trigger
172 if ((chopperTimestamp == 0) && m_discardPreChopperEvents) {
173 // throw away events before first chopper trigger
174 continue;
175 }
176 spectrum.addEventQuickly(Types::Event::TofEvent(double(event.timestamp - chopperTimestamp) / 10.0));
177 }
178
180 oversizedChanelIndexCounterA += oversizedChanelIndexCounter;
181 oversizedPosCounterA += oversizedPosCounter;
182 }
184
185 if (oversizedChanelIndexCounterA > 0) {
186 g_log.warning() << "Bad chanel indices: " << oversizedChanelIndexCounterA << std::endl;
187 }
188 if (oversizedPosCounterA > 0) {
189 g_log.warning() << "Bad position values: " << oversizedPosCounterA << std::endl;
190 }
191}
192
193std::vector<uint8_t> LoadDNSEvent::parse_Header(FileByteStream &file) {
194 // using Boyer-Moore String Search:
195
196 // search for header_sep and store actual header:
197 std::vector<uint8_t> header;
198
199 std::array<uint8_t, header_sep.size()> current_window;
200 file.readRaw(current_window);
201 while (!file.eof()) {
202 if (current_window == header_sep) {
203 return header;
204 } else {
205 const auto orig_header_size = header.size();
206 header.resize(header.size() + 1);
207 const auto win_data = current_window.data();
208 std::copy(win_data, win_data + 1, header.data() + orig_header_size);
209 const std::array<uint8_t, header_sep.size()> orig_window = current_window;
210 file.readRaw(current_window, 1);
211 std::copy(orig_window.data() + 1, orig_window.data() + header_sep.size(), win_data);
212 }
213 }
214 return header;
215}
216
217std::vector<std::vector<uint8_t>> LoadDNSEvent::split_File(FileByteStream &file, const unsigned maxChunckCount) {
218
219 const uint64_t minChunckSize = MAX_BUFFER_BYTES_SIZE;
220 const uint64_t chunckSize = std::max(minChunckSize, file.fileSize() / maxChunckCount);
221
222 std::vector<std::vector<uint8_t>> result;
223
224 while (!file.eof()) {
225 result.push_back(std::vector<uint8_t>());
226 // read a big chunck of file:
227 auto &data = result.back();
228 data.resize(chunckSize);
229 try {
230 file.readRaw(*data.begin(), chunckSize);
231 } catch (std::ifstream::failure &) {
232 data.resize(file.gcount());
233 return result;
234 }
235
236 // search for a block_separator, and append everything up to it :
237 static const auto windowSize = block_sep.size();
238 uint8_t *current_window = nullptr;
239 std::array<uint8_t, windowSize> *windowAsArray =
240 reinterpret_cast<std::array<uint8_t, windowSize> *>(current_window);
241
242 try {
243 data.resize(data.size() + windowSize);
244 // accomodate for possible relocation of vector...:
245 current_window = &*(data.end() - windowSize);
246 windowAsArray = reinterpret_cast<std::array<uint8_t, windowSize> *>(current_window);
247 file.readRaw(current_window[0], windowSize);
248
249 while ((*windowAsArray != block_sep) & (!file.eof())) {
250 const auto orig_data_size = data.size();
251 data.resize(orig_data_size + 1);
252 // accomodate for possible relocation of vector...:
253 current_window = (&data.back() - windowSize + 1);
254 windowAsArray = reinterpret_cast<std::array<uint8_t, windowSize> *>(current_window);
255 file.readRaw(current_window[windowSize - 1], 1);
256 }
257 } catch (std::ifstream::failure &) {
258 return result;
259 }
260 } // while
261 return result;
262}
263
265 // File := Header Body
266 std::vector<uint8_t> header = parse_Header(file);
267
268 // check it is actually a mesytec psd listmode file:
269 if (!startsWith(header, std::string("mesytec psd listmode data"))) {
270 g_log.error() << "This seems not to be a mesytec psd listmode file: " << fileName;
271 throw Exception::FileError("This seems not to be a mesytec psd listmode file: ", fileName);
272 }
273
274 // Parse actual data:
275 const int threadCount = PARALLEL_GET_MAX_THREADS;
276 // Split File:
277 std::vector<std::vector<uint8_t>> filechuncks = split_File(file, threadCount);
278 g_log.debug() << "filechuncks count = " << filechuncks.size() << std::endl;
279
280 std::vector<EventAccumulator> eventAccumulators(filechuncks.size());
281 for (auto &evtAcc : eventAccumulators) {
282 evtAcc.neutronEvents.resize(m_detectorPixelCount);
283 }
284
285 // parse individual file chuncks:
287 for (int i = 0; i < static_cast<int>(filechuncks.size()); ++i) {
288 auto filechunck = filechuncks[static_cast<size_t>(i)];
289 g_log.debug() << "filechunck.size() = " << filechunck.size() << std::endl;
290 auto vbs = VectorByteStream(filechunck);
291 parse_BlockList(vbs, eventAccumulators[static_cast<size_t>(i)]);
292 }
293
294 EventAccumulator finalEventAccumulator;
295 finalEventAccumulator.neutronEvents.resize(m_detectorPixelCount);
296
297 // combine eventAccumulators:
298
299 // combine triggerEvents:
300 for (const auto &v : eventAccumulators) {
301 finalEventAccumulator.triggerEvents.insert(finalEventAccumulator.triggerEvents.end(), v.triggerEvents.begin(),
302 v.triggerEvents.end());
303 }
304
305 // combine neutronEvents:
307 for (int i = 0; i < static_cast<int>(finalEventAccumulator.neutronEvents.size()); ++i) {
308 auto &allNeutronEvents = finalEventAccumulator.neutronEvents[static_cast<size_t>(i)];
309 for (const auto &v : eventAccumulators) {
310 allNeutronEvents.insert(allNeutronEvents.end(), v.neutronEvents[static_cast<size_t>(i)].begin(),
311 v.neutronEvents[static_cast<size_t>(i)].end());
312 }
313 }
314
315 return finalEventAccumulator;
316}
317
319 // BlockList := DataBuffer BlockListTrail
320 while (!file.eof() && file.peek() != 0xFF) {
321 parse_Block(file, eventAccumulator);
322 }
323}
324
326 // Block := DataBufferHeader DataBuffer
327 parse_DataBuffer(file, eventAccumulator);
329}
330
332 auto separator = file.readRaw(separator_t());
333 if (separator != block_sep) {
334 throw std::runtime_error(std::string("File Integrety LOST. 0x"));
335 }
336}
337
339 const auto bufferHeader = parse_DataBufferHeader(file);
340
341 const uint16_t dataLength = uint16_t(bufferHeader.bufferLength - 21);
342 const auto event_count = dataLength / 3;
343
344 for (uint16_t i = 0; i < event_count; i++) {
345 parse_andAddEvent(file, bufferHeader, eventAccumulator);
346 }
347}
348
350 uint16_t ts1 = 0;
351 uint16_t ts2 = 0;
352 uint16_t ts3 = 0;
353 BufferHeader header = {};
354 file.read<2>(header.bufferLength);
355 file.read<2>(header.bufferVersion);
356 file.read<2>(header.headerLength);
357 file.read<2>(header.bufferNumber);
358 file.read<2>(header.runId);
359 file.read<1>(header.mcpdId);
360 file.read<1>(header.deviceStatus);
361 file.read<2>(ts1);
362 file.read<2>(ts2);
363 file.read<2>(ts3);
364 // 48 bit timestamp is 3 2-byte MSB words ordered LSB
365 header.timestamp = (uint64_t)ts3 << 32 | (uint64_t)ts2 << 16 | ts1;
366 file.skip<24>();
367 return header;
368}
369
371 const LoadDNSEvent::BufferHeader &bufferHeader) {
372 TriggerEvent triggerEvent = {};
373 uint8_t trigId = 0;
374 uint8_t dataId = 0;
375 trigId = (data >> 44) & 0b111; // 3 bit
376 dataId = (data >> 40) & 0b1111; // 4 bit
377 triggerEvent.event.timestamp = data & 0x7ffff; // 19bit
378 triggerEvent.event.timestamp += bufferHeader.timestamp;
379 triggerEvent.isChopperTrigger = ((dataId == m_chopperChannel - 1) && (trigId == 7));
380 return triggerEvent;
381}
382
384 const LoadDNSEvent::BufferHeader &bufferHeader) {
385 NeutronEvent neutronEvent = {};
386 uint16_t position = 0;
387 uint8_t modid = 0;
388 uint8_t slotid = 0;
389 modid = (data >> 44) & 0b111; // 3bit
390 slotid = (data >> 39) & 0b11111; // 5bit
391 position = (data >> 19) & 0x3ff; // 10bit
392 neutronEvent.event.timestamp = (data & 0x7ffff); // 19bit
393 neutronEvent.event.timestamp += bufferHeader.timestamp;
394 neutronEvent.wsIndex = (bufferHeader.mcpdId << 6u | modid << 3u | slotid) << 10u | position;
395 return neutronEvent;
396}
397
399 LoadDNSEvent::EventAccumulator &eventAccumulator) {
400 uint16_t data1 = 0;
401 uint16_t data2 = 0;
402 uint16_t data3 = 0;
403 uint64_t data = 0;
404 event_id_e eventId;
405 file.read<2>(data1);
406 file.read<2>(data2);
407 file.read<2>(data3);
408 // 48 bit event is 3 2-byte MSB words ordered LSB
409 data = (uint64_t)data3 << 32 | (uint64_t)data2 << 16 | data1;
410 eventId = static_cast<event_id_e>((data >> 47) & 0b1);
411 switch (eventId) {
412 case event_id_e::TRIGGER: {
413 const auto triggerEvent = processTrigger(data, bufferHeader);
414 if (triggerEvent.isChopperTrigger) {
415 eventAccumulator.triggerEvents.push_back(triggerEvent.event);
416 }
417 } break;
418 case event_id_e::NEUTRON: {
419 const auto neutronEvent = processNeutron(data, bufferHeader);
420 eventAccumulator.neutronEvents[neutronEvent.wsIndex].push_back(neutronEvent.event);
421 } break;
422 default:
423 g_log.error() << "unknown event id " << eventId << "\n";
424 break;
425 }
426}
427
429 auto separator = file.readRaw(separator_t());
430 if (separator != closing_sig) {
431 throw std::runtime_error(std::string("File Integrety LOST. 0x"));
432 }
433}
434
435} // namespace Mantid::DataHandling
double position
Definition: GetAllEi.cpp:154
static constexpr separator_t closing_sig
static constexpr separator_t block_sep
std::array< uint8_t, 8 > separator_t
static constexpr separator_t header_sep
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_FOR_NO_WSP_CHECK()
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_GET_MAX_THREADS
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
bool eof() const
Definition: BitStream.h:61
FileByteStream & readRaw(T &result, const std::size_t &bytecount)
Definition: BitStream.h:41
uint64_t fileSize() const
Definition: BitStream.h:71
std::streamsize gcount() const
Definition: BitStream.h:68
CancelException is thrown to cancel execution of the algorithm.
Definition: Algorithm.h:145
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
bool getCancel() const
Returns the cancellation state.
Definition: Algorithm.cpp:1657
@ 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.
TriggerEvent processTrigger(const uint64_t &data, const LoadDNSEvent::BufferHeader &bufferHeader)
std::vector< uint8_t > parse_Header(FileByteStream &file)
void parse_BlockList(VectorByteStream &file, EventAccumulator &eventAccumulator)
void populate_EventWorkspace(Mantid::DataObjects::EventWorkspace_sptr &eventWS, EventAccumulator &finalEventAccumulator)
void parse_andAddEvent(VectorByteStream &file, const BufferHeader &bufferHeader, EventAccumulator &eventAccumulator)
NeutronEvent processNeutron(const uint64_t &data, const LoadDNSEvent::BufferHeader &bufferHeader)
void exec() override
Run the algorithm.
EventAccumulator parse_File(FileByteStream &file, const std::string &fileName)
void parse_BlockSeparator(VectorByteStream &file)
BufferHeader parse_DataBufferHeader(VectorByteStream &file)
void parse_EndSignature(FileByteStream &file)
void parse_DataBuffer(VectorByteStream &file, EventAccumulator &eventAccumulator)
std::vector< std::vector< uint8_t > > split_File(FileByteStream &file, const unsigned maxChunckCount)
void init() override
Virtual method - must be overridden by concrete algorithm.
void parse_Block(VectorByteStream &file, EventAccumulator &eventAccumulator)
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
Records the filename and the description of failure.
Definition: Exception.h:98
Defines a wrapper around an open file.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
VectorByteStream & read(T &result)
Definition: BitStream.h:117
uint8_t peek()
Definition: BitStream.h:123
VectorByteStream & skip()
Definition: BitStream.h:127
VectorByteStream & readRaw(T &result, const std::size_t &bytecount)
Definition: BitStream.h:98
bool eof() const
Definition: BitStream.h:125
Kernel::Logger g_log("ExperimentInfo")
static logger object
bool startsWith(const string &str, const string &prefix)
Returns true if str starts with prefix.
const unsigned int PIXEL_PER_TUBE
const unsigned int MAX_BUFFER_BYTES_SIZE
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
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.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
std::vector< std::vector< CompactEvent > > neutronEvents
Definition: LoadDNSEvent.h:77
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54