ProcessBankData.cpp 11.2 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
David Fairbrother's avatar
David Fairbrother committed
7
8
#include <utility>

9
#include "MantidDataHandling/DefaultEventLoader.h"
10
#include "MantidDataHandling/LoadEventNexus.h"
David Fairbrother's avatar
David Fairbrother committed
11
#include "MantidDataHandling/ProcessBankData.h"
Alex Buts's avatar
Alex Buts committed
12

13
using namespace Mantid::DataObjects;
14

15
namespace Mantid::DataHandling {
16

Samuel Jones's avatar
Samuel Jones committed
17
18
19
20
21
22
23
ProcessBankData::ProcessBankData(DefaultEventLoader &m_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)
David Fairbrother's avatar
David Fairbrother committed
24
    : Task(), m_loader(m_loader), entry_name(std::move(entry_name)),
Samuel Jones's avatar
Samuel Jones committed
25
26
27
28
29
      pixelID_to_wi_vector(m_loader.pixelID_to_wi_vector), pixelID_to_wi_offset(m_loader.pixelID_to_wi_offset),
      prog(prog), event_id(std::move(event_id)), event_time_of_flight(std::move(event_time_of_flight)),
      numEvents(numEvents), startAt(startAt), event_index(std::move(event_index)),
      thisBankPulseTimes(std::move(thisBankPulseTimes)), have_weight(have_weight),
      event_weight(std::move(event_weight)), m_min_id(min_event_id), m_max_id(max_event_id) {
30
31
32
  // Cost is approximately proportional to the number of events to process.
  m_cost = static_cast<double>(numEvents);
}
33

34
35
36
namespace {
// this assumes that last_pulse_index is already to the point of including this
// one so we only need to search forward
Samuel Jones's avatar
Samuel Jones committed
37
38
inline size_t getPulseIndex(const size_t event_index, const size_t last_pulse_index,
                            const std::shared_ptr<std::vector<uint64_t>> &event_index_vec) {
Peterson, Peter's avatar
Peterson, Peter committed
39
  if (last_pulse_index + 1 >= event_index_vec->size())
40
41
42
43
44
45
46
47
    return last_pulse_index;

  // linear search is used because it is more likely that the next pulse index
  // is the correct one to use the + last_pulse_index + 1 is because we are
  // confirm that the next index is bigger, not the current
  const auto event_index_end = event_index_vec->cend();
  auto event_index_iter = event_index_vec->cbegin() + last_pulse_index;

Samuel Jones's avatar
Samuel Jones committed
48
  while ((event_index < *event_index_iter) || (event_index >= *(event_index_iter + 1))) {
49
50
51
52
53
54
55
56
57
58
    event_index_iter++;

    // make sure not to go past the end
    if (event_index_iter + 1 == event_index_end)
      break;
  }
  return std::distance(event_index_vec->cbegin(), event_index_iter);
}
} // namespace

59
60
/** Run the data processing
 * FIXME/TODO - split run() into readable methods
61
 */
62
void ProcessBankData::run() { // override {
63
  // Local tof limits
Samuel Jones's avatar
Samuel Jones committed
64
  double my_shortest_tof = static_cast<double>(std::numeric_limits<uint32_t>::max()) * 0.1;
65
66
67
68
  double my_longest_tof = 0.;
  // A count of "bad" TOFs that were too high
  size_t badTofs = 0;
  size_t my_discarded_events(0);
69

70
71
  prog->report(entry_name + ": precount");
  // ---- Pre-counting events per pixel ID ----
72
73
74
  auto &outputWS = m_loader.m_ws;
  auto *alg = m_loader.alg;
  if (m_loader.precount) {
75

76
77
    std::vector<size_t> counts(m_max_id - m_min_id + 1, 0);
    for (size_t i = 0; i < numEvents; i++) {
78
      const auto thisId = detid_t((*event_id)[i]);
79
80
81
      if (thisId >= m_min_id && thisId <= m_max_id)
        counts[thisId - m_min_id]++;
    }
82

83
84
85
    // Now we pre-allocate (reserve) the vectors of events in each pixel
    // counted
    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
90
91
92
        // Find the the workspace index corresponding to that pixel ID
        // Allocate it
        if (wi < numEventLists) {
          outputWS.reserveEventListAt(wi, counts[pixID - m_min_id]);
93
        }
94
95
        if (alg->getCancel())
          break; // User cancellation
96
97
      }
    }
98
  }
99

100
101
102
103
  // Check for canceled algorithm
  if (alg->getCancel()) {
    return;
  }
104

105
  // Default pulse time (if none are found)
Samuel Jones's avatar
Samuel Jones committed
106
  const bool pulsetimesincreasing =
107
      std::is_sorted(thisBankPulseTimes->pulseTimes.cbegin(), thisBankPulseTimes->pulseTimes.cend());
Peterson, Peter's avatar
Peterson, Peter committed
108
109
  if (!std::is_sorted(event_index->cbegin(), event_index->cend()))
    throw std::runtime_error("Event index is not sorted");
110

111
  // And there are this many pulses
112
  const auto NUM_PULSES = thisBankPulseTimes->pulseTimes.size();
113
  prog->report(entry_name + ": filling events");
114

115
  // Will we need to compress?
116
  const bool compress = (alg->compressTolerance >= 0);
117

118
  // Which detector IDs were touched?
119
  std::vector<bool> usedDetIds;
120
  usedDetIds.assign(m_max_id - m_min_id + 1, false);
121

Peterson, Peter's avatar
Peterson, Peter committed
122
123
  const double TOF_MIN = alg->filter_tof_min;
  const double TOF_MAX = alg->filter_tof_max;
Peterson, Peter's avatar
Peterson, Peter committed
124

Samuel Jones's avatar
Samuel Jones committed
125
  for (std::size_t pulseIndex = getPulseIndex(startAt, 0, event_index); pulseIndex < NUM_PULSES; pulseIndex++) {
Peterson, Peter's avatar
Peterson, Peter committed
126
127
128
129
    // Save the pulse time at this index for creating those events
    const auto pulsetime = thisBankPulseTimes->pulseTimes[pulseIndex];
    const int logPeriodNumber = thisBankPulseTimes->periodNumbers[pulseIndex];
    const int periodIndex = logPeriodNumber - 1;
130

Peterson, Peter's avatar
Peterson, Peter committed
131
132
133
134
135
136
137
138
139
    const auto firstEventIndex = getFirstEventIndex(pulseIndex);
    if (firstEventIndex > numEvents)
      break;

    const auto lastEventIndex = getLastEventIndex(pulseIndex, NUM_PULSES);
    if (firstEventIndex == lastEventIndex)
      continue;
    else if (firstEventIndex > lastEventIndex) {
      std::stringstream msg;
Samuel Jones's avatar
Samuel Jones committed
140
141
      msg << "Something went really wrong: " << firstEventIndex << " > " << lastEventIndex << "| " << entry_name
          << " startAt=" << startAt << " numEvents=" << event_index->size() << " RAWINDICES=["
Peterson, Peter's avatar
Peterson, Peter committed
142
143
144
          << firstEventIndex + startAt << ",?)"
          << " pulseIndex=" << pulseIndex << " of " << event_index->size();
      throw std::runtime_error(msg.str());
Peterson, Peter's avatar
Peterson, Peter committed
145
    }
146

Samuel Jones's avatar
Samuel Jones committed
147
    for (std::size_t eventIndex = firstEventIndex; eventIndex < lastEventIndex; ++eventIndex) {
Peterson, Peter's avatar
Peterson, Peter committed
148
149
      // We cached a pointer to the vector<tofEvent> -> so retrieve it and add
      // the event
150
      const detid_t detId = (*event_id)[eventIndex];
Peterson, Peter's avatar
Peterson, Peter committed
151
152
      if (detId >= m_min_id && detId <= m_max_id) {
        // Create the tofevent
Samuel Jones's avatar
Samuel Jones committed
153
        const auto tof = static_cast<double>((*event_time_of_flight)[eventIndex]);
Peterson, Peter's avatar
Peterson, Peter committed
154
155
156
157
        // this is fancy for check if value is in range
        if ((tof - TOF_MIN) * (tof - TOF_MAX) <= 0.) {
          // Handle simulated data if present
          if (have_weight) {
Samuel Jones's avatar
Samuel Jones committed
158
            auto *eventVector = m_loader.weightedEventVectors[periodIndex][detId];
Peterson, Peter's avatar
Peterson, Peter committed
159
160
            // NULL eventVector indicates a bad spectrum lookup
            if (eventVector) {
Samuel Jones's avatar
Samuel Jones committed
161
              const auto weight = static_cast<double>((*event_weight)[eventIndex]);
Peterson, Peter's avatar
Peterson, Peter committed
162
163
164
165
166
              const double errorSq = weight * weight;
              eventVector->emplace_back(tof, pulsetime, weight, errorSq);
            } else {
              ++my_discarded_events;
            }
167
          } else {
Peterson, Peter's avatar
Peterson, Peter committed
168
169
170
171
172
173
174
175
            // We have cached the vector of events for this detector ID
            auto *eventVector = m_loader.eventVectors[periodIndex][detId];
            // NULL eventVector indicates a bad spectrum lookup
            if (eventVector) {
              eventVector->emplace_back(tof, pulsetime);
            } else {
              ++my_discarded_events;
            }
176
          }
177

Peterson, Peter's avatar
Peterson, Peter committed
178
179
180
181
182
183
184
185
186
187
188
189
          // Skip any events that are the cause of bad DAS data (e.g. a negative
          // number in uint32 -> 2.4 billion * 100 nanosec = 2.4e8 microsec)
          if (tof < 2e8) {
            // tof limits from things observed here
            if (tof > my_longest_tof) {
              my_longest_tof = tof;
            }
            if (tof < my_shortest_tof) {
              my_shortest_tof = tof;
            }
          } else
            badTofs++;
190

191
192
          // Track all the touched wi
          usedDetIds[detId - m_min_id] = true;
Peterson, Peter's avatar
Peterson, Peter committed
193
        } // valid time-of-flight
194

Peterson, Peter's avatar
Peterson, Peter committed
195
196
197
198
199
200
201
202
203
204
205
      } // valid detector IDs
    }   // for events in pulse
    // check if cancelled after each pulse
    if (alg->getCancel())
      break;
  } // for pulses

  // Check for canceled algorithm
  if (alg->getCancel()) {
    return;
  }
206

207
208
  //------------ Compress Events (or set sort order) ------------------
  // Do it on all the detector IDs we touched
209
  const size_t numEventLists = outputWS.getNumberHistograms();
210
211
212
213
  for (detid_t pixID = m_min_id; pixID <= m_max_id; ++pixID) {
    if (usedDetIds[pixID - m_min_id]) {
      // Find the the workspace index corresponding to that pixel ID
      size_t wi = getWorkspaceIndexFromPixelID(pixID);
214
215
216
217
218
219
220
221
222
223
      if (wi < numEventLists) {
        auto &el = outputWS.getSpectrum(wi);
        if (compress)
          el.compressEvents(alg->compressTolerance, &el);
        else {
          if (pulsetimesincreasing)
            el.setSortOrder(DataObjects::PULSETIME_SORT);
          else
            el.setSortOrder(DataObjects::UNSORTED);
        }
224
225
      }
    }
226
227
  }
  prog->report(entry_name + ": filled events");
228

Samuel Jones's avatar
Samuel Jones committed
229
  alg->getLogger().debug() << entry_name << (pulsetimesincreasing ? " had " : " DID NOT have ")
230
                           << "monotonically increasing pulse times\n";
231

232
233
234
235
236
237
238
239
240
  // Join back up the tof limits to the global ones
  // This is not thread safe, so only one thread at a time runs this.
  {
    std::lock_guard<std::mutex> _lock(alg->m_tofMutex);
    if (my_shortest_tof < alg->shortest_tof) {
      alg->shortest_tof = my_shortest_tof;
    }
    if (my_longest_tof > alg->longest_tof) {
      alg->longest_tof = my_longest_tof;
241
    }
242
243
244
    alg->bad_tofs += badTofs;
    alg->discarded_events += my_discarded_events;
  }
245
246

#ifndef _WIN32
Samuel Jones's avatar
Samuel Jones committed
247
  alg->getLogger().debug() << "Time to process " << entry_name << " " << m_timer << "\n";
248
#endif
249
} // END-OF-RUN()
250

Peterson, Peter's avatar
Peterson, Peter committed
251
252
253
254
255
256
257
258
size_t ProcessBankData::getFirstEventIndex(const size_t pulseIndex) const {
  const auto firstEventIndex = event_index->operator[](pulseIndex);
  if (firstEventIndex >= startAt)
    return firstEventIndex - startAt;
  else
    return 0;
}

Samuel Jones's avatar
Samuel Jones committed
259
size_t ProcessBankData::getLastEventIndex(const size_t pulseIndex, const size_t numPulses) const {
Peterson, Peter's avatar
Peterson, Peter committed
260
261
262
  if (pulseIndex + 1 >= numPulses)
    return numEvents;

Samuel Jones's avatar
Samuel Jones committed
263
  const size_t lastEventIndex = event_index->operator[](pulseIndex + 1) - startAt;
Peterson, Peter's avatar
Peterson, Peter committed
264
265
266
267

  return std::min(lastEventIndex, numEvents);
}

268
269
270
271
/**
 * Get the workspace index for a given pixel ID. Throws if the pixel ID is
 * not in the expected range.
 *
272
 * @param pixID :: The pixel ID to look up
273
274
 * @return The workspace index for this pixel
 */
275
276
277
size_t ProcessBankData::getWorkspaceIndexFromPixelID(const detid_t pixID) {
  // Check that the vector index is not out of range
  const detid_t offset_pixID = pixID + pixelID_to_wi_offset;
Samuel Jones's avatar
Samuel Jones committed
278
  if (offset_pixID < 0 || offset_pixID >= static_cast<int32_t>(pixelID_to_wi_vector.size())) {
279
    std::stringstream msg;
Samuel Jones's avatar
Samuel Jones committed
280
    msg << "Error finding workspace index; pixelID " << pixID << " with offset " << pixelID_to_wi_offset
Gemma Guest's avatar
Gemma Guest committed
281
        << " is out of range (length=" << pixelID_to_wi_vector.size() << ")";
282
283
284
285
    throw std::runtime_error(msg.str());
  }
  return pixelID_to_wi_vector[offset_pixID];
}
286
} // namespace Mantid::DataHandling