Mantid
Loading...
Searching...
No Matches
ProcessBankData.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 +
7#include <utility>
8
12
13using namespace Mantid::DataObjects;
14
15namespace Mantid::DataHandling {
16
17ProcessBankData::ProcessBankData(DefaultEventLoader &m_loader, std::string entry_name, API::Progress *prog,
18 std::shared_ptr<std::vector<uint32_t>> event_id,
19 std::shared_ptr<std::vector<float>> event_time_of_flight, size_t numEvents,
20 size_t startAt, std::shared_ptr<std::vector<uint64_t>> event_index,
21 std::shared_ptr<BankPulseTimes> thisBankPulseTimes, bool have_weight,
22 std::shared_ptr<std::vector<float>> event_weight, detid_t min_event_id,
23 detid_t max_event_id)
24 : Task(), m_loader(m_loader), entry_name(std::move(entry_name)),
25 pixelID_to_wi_vector(m_loader.pixelID_to_wi_vector), pixelID_to_wi_offset(m_loader.pixelID_to_wi_offset),
26 prog(prog), event_id(std::move(event_id)), event_time_of_flight(std::move(event_time_of_flight)),
27 numEvents(numEvents), startAt(startAt), event_index(std::move(event_index)),
28 thisBankPulseTimes(std::move(thisBankPulseTimes)), have_weight(have_weight),
29 event_weight(std::move(event_weight)), m_min_id(min_event_id), m_max_id(max_event_id) {
30 // Cost is approximately proportional to the number of events to process.
31 m_cost = static_cast<double>(numEvents);
32}
33
34namespace {
35// this assumes that last_pulse_index is already to the point of including this
36// one so we only need to search forward
37inline size_t getPulseIndex(const size_t event_index, const size_t last_pulse_index,
38 const std::shared_ptr<std::vector<uint64_t>> &event_index_vec) {
39 if (last_pulse_index + 1 >= event_index_vec->size())
40 return last_pulse_index;
41
42 // linear search is used because it is more likely that the next pulse index
43 // is the correct one to use the + last_pulse_index + 1 is because we are
44 // confirm that the next index is bigger, not the current
45 const auto event_index_end = event_index_vec->cend();
46 auto event_index_iter = event_index_vec->cbegin() + last_pulse_index;
47
48 while ((event_index < *event_index_iter) || (event_index >= *(event_index_iter + 1))) {
49 event_index_iter++;
50
51 // make sure not to go past the end
52 if (event_index_iter + 1 == event_index_end)
53 break;
54 }
55 return std::distance(event_index_vec->cbegin(), event_index_iter);
56}
57} // namespace
58
62void ProcessBankData::run() { // override {
63 // Local tof limits
64 double my_shortest_tof = static_cast<double>(std::numeric_limits<uint32_t>::max()) * 0.1;
65 double my_longest_tof = 0.;
66 // A count of "bad" TOFs that were too high
67 size_t badTofs = 0;
68 size_t my_discarded_events(0);
69
70 prog->report(entry_name + ": precount");
71 // ---- Pre-counting events per pixel ID ----
72 auto &outputWS = m_loader.m_ws;
73 auto *alg = m_loader.alg;
74 if (m_loader.precount) {
75
76 std::vector<size_t> counts(m_max_id - m_min_id + 1, 0);
77 for (size_t i = 0; i < numEvents; i++) {
78 const auto thisId = detid_t((*event_id)[i]);
79 if (thisId >= m_min_id && thisId <= m_max_id)
80 counts[thisId - m_min_id]++;
81 }
82
83 // Now we pre-allocate (reserve) the vectors of events in each pixel
84 // counted
85 const size_t numEventLists = outputWS.getNumberHistograms();
86 for (detid_t pixID = m_min_id; pixID <= m_max_id; ++pixID) {
87 if (counts[pixID - m_min_id] > 0) {
88 size_t wi = getWorkspaceIndexFromPixelID(pixID);
89 // Find the the workspace index corresponding to that pixel ID
90 // Allocate it
91 if (wi < numEventLists) {
92 outputWS.reserveEventListAt(wi, counts[pixID - m_min_id]);
93 }
94 if (alg->getCancel())
95 return; // User cancellation
96 }
97 }
98 }
99
100 // Default pulse time (if none are found)
101 const bool pulsetimesincreasing =
102 std::is_sorted(thisBankPulseTimes->pulseTimes.cbegin(), thisBankPulseTimes->pulseTimes.cend());
103 if (!std::is_sorted(event_index->cbegin(), event_index->cend()))
104 throw std::runtime_error("Event index is not sorted");
105
106 // And there are this many pulses
107 const auto NUM_PULSES = thisBankPulseTimes->pulseTimes.size();
108 prog->report(entry_name + ": filling events");
109
110 // Will we need to compress?
111 const bool compress = (alg->compressTolerance >= 0);
112
113 // Which detector IDs were touched?
114 std::vector<bool> usedDetIds;
115 usedDetIds.assign(m_max_id - m_min_id + 1, false);
116
117 const double TOF_MIN = alg->filter_tof_min;
118 const double TOF_MAX = alg->filter_tof_max;
119
120 for (std::size_t pulseIndex = getPulseIndex(startAt, 0, event_index); pulseIndex < NUM_PULSES; pulseIndex++) {
121 // Save the pulse time at this index for creating those events
122 const auto pulsetime = thisBankPulseTimes->pulseTimes[pulseIndex];
123 const int logPeriodNumber = thisBankPulseTimes->periodNumbers[pulseIndex];
124 const int periodIndex = logPeriodNumber - 1;
125
126 const auto firstEventIndex = getFirstEventIndex(pulseIndex);
127 if (firstEventIndex > numEvents)
128 break;
129
130 const auto lastEventIndex = getLastEventIndex(pulseIndex, NUM_PULSES);
131 if (firstEventIndex == lastEventIndex)
132 continue;
133 else if (firstEventIndex > lastEventIndex) {
134 std::stringstream msg;
135 msg << "Something went really wrong: " << firstEventIndex << " > " << lastEventIndex << "| " << entry_name
136 << " startAt=" << startAt << " numEvents=" << event_index->size() << " RAWINDICES=["
137 << firstEventIndex + startAt << ",?)"
138 << " pulseIndex=" << pulseIndex << " of " << event_index->size();
139 throw std::runtime_error(msg.str());
140 }
141
142 for (std::size_t eventIndex = firstEventIndex; eventIndex < lastEventIndex; ++eventIndex) {
143 // We cached a pointer to the vector<tofEvent> -> so retrieve it and add
144 // the event
145 const detid_t detId = (*event_id)[eventIndex];
146 if (detId >= m_min_id && detId <= m_max_id) {
147 // Create the tofevent
148 const auto tof = static_cast<double>((*event_time_of_flight)[eventIndex]);
149 // this is fancy for check if value is in range
150 if ((tof - TOF_MIN) * (tof - TOF_MAX) <= 0.) {
151 // Handle simulated data if present
152 if (have_weight) {
153 auto *eventVector = m_loader.weightedEventVectors[periodIndex][detId];
154 // NULL eventVector indicates a bad spectrum lookup
155 if (eventVector) {
156 const auto weight = static_cast<double>((*event_weight)[eventIndex]);
157 const double errorSq = weight * weight;
158 eventVector->emplace_back(tof, pulsetime, weight, errorSq);
159 } else {
160 ++my_discarded_events;
161 }
162 } else {
163 // We have cached the vector of events for this detector ID
164 auto *eventVector = m_loader.eventVectors[periodIndex][detId];
165 // NULL eventVector indicates a bad spectrum lookup
166 if (eventVector) {
167 eventVector->emplace_back(tof, pulsetime);
168 } else {
169 ++my_discarded_events;
170 }
171 }
172
173 // Skip any events that are the cause of bad DAS data (e.g. a negative
174 // number in uint32 -> 2.4 billion * 100 nanosec = 2.4e8 microsec)
175 if (tof < 2e8) {
176 // tof limits from things observed here
177 if (tof > my_longest_tof) {
178 my_longest_tof = tof;
179 }
180 if (tof < my_shortest_tof) {
181 my_shortest_tof = tof;
182 }
183 } else
184 badTofs++;
185
186 // Track all the touched wi
187 usedDetIds[detId - m_min_id] = true;
188 } // valid time-of-flight
189
190 } // valid detector IDs
191 } // for events in pulse
192 // check if cancelled after each pulse
193 if (alg->getCancel())
194 return;
195 } // for pulses
196
197 //------------ Compress Events (or set sort order) ------------------
198 // Do it on all the detector IDs we touched
199 const size_t numEventLists = outputWS.getNumberHistograms();
200 for (detid_t pixID = m_min_id; pixID <= m_max_id; ++pixID) {
201 if (usedDetIds[pixID - m_min_id]) {
202 // Find the the workspace index corresponding to that pixel ID
203 size_t wi = getWorkspaceIndexFromPixelID(pixID);
204 if (wi < numEventLists) {
205 auto &el = outputWS.getSpectrum(wi);
206 if (compress)
207 el.compressEvents(alg->compressTolerance, &el);
208 else {
209 if (pulsetimesincreasing)
210 el.setSortOrder(DataObjects::PULSETIME_SORT);
211 else
212 el.setSortOrder(DataObjects::UNSORTED);
213 }
214 }
215 }
216 }
217 prog->report(entry_name + ": filled events");
218
219 alg->getLogger().debug() << entry_name << (pulsetimesincreasing ? " had " : " DID NOT have ")
220 << "monotonically increasing pulse times\n";
221
222 // Join back up the tof limits to the global ones
223 // This is not thread safe, so only one thread at a time runs this.
224 {
225 std::lock_guard<std::mutex> _lock(alg->m_tofMutex);
226 if (my_shortest_tof < alg->shortest_tof) {
227 alg->shortest_tof = my_shortest_tof;
228 }
229 if (my_longest_tof > alg->longest_tof) {
230 alg->longest_tof = my_longest_tof;
231 }
232 alg->bad_tofs += badTofs;
233 alg->discarded_events += my_discarded_events;
234 }
235
236#ifndef _WIN32
237 alg->getLogger().debug() << "Time to process " << entry_name << " " << m_timer << "\n";
238#endif
239} // END-OF-RUN()
240
241size_t ProcessBankData::getFirstEventIndex(const size_t pulseIndex) const {
242 const auto firstEventIndex = event_index->operator[](pulseIndex);
243 if (firstEventIndex >= startAt)
244 return firstEventIndex - startAt;
245 else
246 return 0;
247}
248
249size_t ProcessBankData::getLastEventIndex(const size_t pulseIndex, const size_t numPulses) const {
250 if (pulseIndex + 1 >= numPulses)
251 return numEvents;
252
253 const size_t lastEventIndex = event_index->operator[](pulseIndex + 1) - startAt;
254
255 return std::min(lastEventIndex, numEvents);
256}
257
266 // Check that the vector index is not out of range
267 const detid_t offset_pixID = pixID + pixelID_to_wi_offset;
268 if (offset_pixID < 0 || offset_pixID >= static_cast<int32_t>(pixelID_to_wi_vector.size())) {
269 std::stringstream msg;
270 msg << "Error finding workspace index; pixelID " << pixID << " with offset " << pixelID_to_wi_offset
271 << " is out of range (length=" << pixelID_to_wi_vector.size() << ")";
272 throw std::runtime_error(msg.str());
273 }
274 return pixelID_to_wi_vector[offset_pixID];
275}
276} // namespace Mantid::DataHandling
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
Helper class for LoadEventNexus that is specific to the current default loading code for NXevent_data...
std::vector< std::vector< std::vector< Mantid::Types::Event::TofEvent > * > > eventVectors
Vector where index = event_id; value = ptr to std::vector<TofEvent> in the event list.
bool precount
Do we pre-count the # of events in each pixel ID?
std::vector< std::vector< std::vector< Mantid::DataObjects::WeightedEvent > * > > weightedEventVectors
Vector where index = event_id; value = ptr to std::vector<WeightedEvent> in the event list.
size_t getWorkspaceIndexFromPixelID(const detid_t pixID)
Get the workspace index for a given pixel ID.
bool have_weight
Flag for simulated data.
size_t getFirstEventIndex(const size_t pulseIndex) const
Mantid::Kernel::Timer m_timer
timer for performance
DefaultEventLoader & m_loader
Algorithm being run.
detid_t m_max_id
Maximum pixel id.
std::shared_ptr< std::vector< uint64_t > > event_index
vector of event index (length of # of pulses)
std::string entry_name
NXS path to bank.
size_t startAt
index of the first event from event_index
std::shared_ptr< std::vector< uint32_t > > event_id
event pixel ID array
detid_t m_min_id
Minimum pixel id.
detid_t pixelID_to_wi_offset
Offset in the pixelID_to_wi_vector to use.
const std::vector< size_t > & pixelID_to_wi_vector
Vector where (index = pixel ID+pixelID_to_wi_offset), value = workspace index)
API::Progress * prog
Progress reporting.
ProcessBankData(DefaultEventLoader &loader, std::string entry_name, API::Progress *prog, std::shared_ptr< std::vector< uint32_t > > event_id, std::shared_ptr< std::vector< float > > event_time_of_flight, size_t numEvents, size_t startAt, std::shared_ptr< std::vector< uint64_t > > event_index, std::shared_ptr< BankPulseTimes > thisBankPulseTimes, bool have_weight, std::shared_ptr< std::vector< float > > event_weight, detid_t min_event_id, detid_t max_event_id)
Constructor.
std::shared_ptr< BankPulseTimes > thisBankPulseTimes
Pulse times for this bank.
size_t getLastEventIndex(const size_t pulseIndex, const size_t numPulses) const
void run() override
Run the data processing FIXME/TODO - split run() into readable methods.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
A Task is a unit of work to be scheduled and run by a ThreadPool.
Definition: Task.h:29
double m_cost
Cached computational cost for the thread.
Definition: Task.h:82
std::size_t numEvents(::NeXus::File &file, bool &hasTotalCounts, bool &oldNeXusFileNames, const std::string &prefix, const NexusHDF5Descriptor &descriptor)
Get the number of events in the currently opened group.
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
STL namespace.