LoadEventPreNexus.cpp 34 KB
Newer Older
1
#include "MantidDataHandling/LoadEventPreNexus.h"
2
#include "MantidAPI/Axis.h"
3
#include "MantidAPI/FileFinder.h"
4
#include "MantidAPI/RegisterFileLoader.h"
Nick Draper's avatar
re #100    
Nick Draper committed
5
6
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/EventList.h"
7
#include "MantidKernel/ArrayProperty.h"
8
#include "MantidKernel/FileValidator.h"
9
#include "MantidKernel/DateAndTime.h"
10
#include "MantidKernel/Glob.h"
11
#include "MantidAPI/FileProperty.h"
12
#include "MantidKernel/ConfigService.h"
13
#include "MantidKernel/BinaryFile.h"
Gigg, Martyn Anthony's avatar
Gigg, Martyn Anthony committed
14
#include "MantidKernel/InstrumentInfo.h"
15
#include "MantidKernel/System.h"
16
#include "MantidKernel/TimeSeriesProperty.h"
17
#include "MantidKernel/UnitFactory.h"
18
#include "MantidKernel/DateAndTime.h"
19
#include "MantidGeometry/IDetector.h"
20
#include "MantidGeometry/Instrument.h"
21
#include "MantidKernel/CPUTimer.h"
22
#include "MantidKernel/VisibleWhenProperty.h"
23
24
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
25
#include "MantidAPI/MemoryManager.h"
26

27
28
29
30
31
32
33
34
35
36
#include <algorithm>
#include <sstream>
#include <stdexcept>
#include <functional>
#include <set>
#include <vector>
#include <Poco/File.h>
#include <Poco/Path.h>
#include <boost/timer.hpp>

37
38
namespace Mantid {
namespace DataHandling {
39

40
DECLARE_FILELOADER_ALGORITHM(LoadEventPreNexus)
41

Nick Draper's avatar
re #100    
Nick Draper committed
42
43
44
using namespace Kernel;
using namespace API;
using namespace Geometry;
45
using boost::posix_time::ptime;
46
using boost::posix_time::time_duration;
Nick Draper's avatar
re #100    
Nick Draper committed
47
48
49
50
51
52
53
54
55
using DataObjects::EventList;
using DataObjects::EventWorkspace;
using DataObjects::EventWorkspace_sptr;
using DataObjects::TofEvent;
using std::cout;
using std::endl;
using std::ifstream;
using std::runtime_error;
using std::stringstream;
56
57
58
59
using std::string;
using std::vector;

// constants for locating the parameters to use in execution
Nick Draper's avatar
re #100    
Nick Draper committed
60
61
62
static const string EVENT_PARAM("EventFilename");
static const string PULSEID_PARAM("PulseidFilename");
static const string MAP_PARAM("MappingFilename");
63
static const string PID_PARAM("SpectrumList");
64
65
static const string PARALLEL_PARAM("UseParallelProcessing");
static const string BLOCK_SIZE_PARAM("LoadingBlockSize");
66
67
static const string OUT_PARAM("OutputWorkspace");

68
69
70
static const string PULSE_EXT("pulseid.dat");
static const string EVENT_EXT("event.dat");

71
/// All pixel ids with matching this mask are errors.
Nick Draper's avatar
re #100    
Nick Draper committed
72
static const PixelType ERROR_PID = 0x80000000;
73
74
/// The maximum possible tof as native type
static const uint32_t MAX_TOF_UINT32 = std::numeric_limits<uint32_t>::max();
Nick Draper's avatar
re #100    
Nick Draper committed
75

76
77
78
79
/// Conversion factor between 100 nanoseconds and 1 microsecond.
static const double TOF_CONVERSION = .1;
/// Conversion factor between picoColumbs and microAmp*hours
static const double CURRENT_CONVERSION = 1.e-6 / 3600.;
Nick Draper's avatar
re #100    
Nick Draper committed
80

81
LoadEventPreNexus::LoadEventPreNexus()
82
    : Mantid::API::IFileLoader<Kernel::FileDescriptor>(), prog(nullptr),
83
84
      spectra_list(), pulsetimes(), event_indices(), proton_charge(),
      proton_charge_tot(0), pixel_to_wkspindex(), pixelmap(), detid_max(),
85
      eventfile(nullptr), num_events(0), num_pulses(0), numpixel(0),
86
87
88
      num_good_events(0), num_error_events(0), num_ignored_events(0),
      first_event(0), max_events(0), using_mapping_file(false),
      loadOnlySomeSpectra(false), spectraLoadMap(), longest_tof(0),
Nick Draper's avatar
Nick Draper committed
89
      shortest_tof(0), parallelProcessing(false) {
90
91
  this->useAlgorithm("LoadEventPreNexus", 2);
}
92

93
LoadEventPreNexus::~LoadEventPreNexus() { delete this->eventfile; }
94

95
/**
96
97
 * Return the confidence with with this algorithm can load the file
 * @param descriptor A descriptor for the file
98
99
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
100
 */
101
102
103
int LoadEventPreNexus::confidence(Kernel::FileDescriptor &descriptor) const {
  if (descriptor.extension().rfind("dat") == std::string::npos)
    return 0;
104
105
106

  // If this looks like a binary file where the exact file length is a multiple
  // of the DasEvent struct then we're probably okay.
107
108
  if (descriptor.isAscii())
    return 0;
109

110
111
  const size_t objSize = sizeof(DasEvent);
  auto &handle = descriptor.data();
112
113
  // get the size of the file in bytes and reset the handle back to the
  // beginning
114
115
116
117
  handle.seekg(0, std::ios::end);
  const size_t filesize = static_cast<size_t>(handle.tellg());
  handle.seekg(0, std::ios::beg);

118
119
120
121
  if (filesize % objSize == 0)
    return 60;
  else
    return 0;
122
123
}

124
//-----------------------------------------------------------------------------
125
/** Initialize the algorithm */
126
void LoadEventPreNexus::init() {
127
  // which files to use
128
129
130
131
132
133
134
135
136
137
138
139
140
141
  declareProperty(
      new FileProperty(EVENT_PARAM, "", FileProperty::Load, EVENT_EXT),
      "The name of the neutron event file to read, including its full or "
      "relative path. The file typically ends in neutron_event.dat (N.B. case "
      "sensitive if running on Linux).");
  declareProperty(new FileProperty(PULSEID_PARAM, "",
                                   FileProperty::OptionalLoad, PULSE_EXT),
                  "File containing the accelerator pulse information; the "
                  "filename will be found automatically if not specified.");
  declareProperty(
      new FileProperty(MAP_PARAM, "", FileProperty::OptionalLoad, ".dat"),
      "File containing the pixel mapping (DAS pixels to pixel IDs) file "
      "(typically INSTRUMENT_TS_YYYY_MM_DD.dat). The filename will be found "
      "automatically if not specified.");
142
143

  // which pixels to load
144
  declareProperty(new ArrayProperty<int64_t>(PID_PARAM),
145
146
                  "A list of individual spectra (pixel IDs) to read, specified "
                  "as e.g. 10:20. Only used if set.");
147

148
  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
149
150
  mustBePositive->setLower(1);
  declareProperty("ChunkNumber", EMPTY_INT(), mustBePositive,
151
152
                  "If loading the file by sections ('chunks'), this is the "
                  "section number of this execution of the algorithm.");
153
  declareProperty("TotalChunks", EMPTY_INT(), mustBePositive,
154
155
                  "If loading the file by sections ('chunks'), this is the "
                  "total number of sections.");
156
  // TotalChunks is only meaningful if ChunkNumber is set
157
158
159
160
  // Would be nice to be able to restrict ChunkNumber to be <= TotalChunks at
  // validation
  setPropertySettings("TotalChunks",
                      new VisibleWhenProperty("ChunkNumber", IS_NOT_DEFAULT));
161

162
  std::vector<std::string> propOptions{"Auto", "Serial", "Parallel"};
163
164
165
166
167
168
169
  declareProperty("UseParallelProcessing", "Auto",
                  boost::make_shared<StringListValidator>(propOptions),
                  "Use multiple cores for loading the data?\n"
                  "  Auto: Use serial loading for small data sets, parallel "
                  "for large data sets.\n"
                  "  Serial: Use a single core.\n"
                  "  Parallel: Use all available cores.");
170
171

  // the output workspace name
172
173
174
175
  declareProperty(
      new WorkspaceProperty<IEventWorkspace>(OUT_PARAM, "", Direction::Output),
      "The name of the workspace that will be created, filled with the read-in "
      "data and stored in the [[Analysis Data Service]].");
176
177
}

178
//-----------------------------------------------------------------------------
179
static string generatePulseidName(string eventfile) {
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
  size_t start;
  string ending;

  // normal ending
  ending = "neutron_event.dat";
  start = eventfile.find(ending);
  if (start != string::npos)
    return eventfile.replace(start, ending.size(), "pulseid.dat");

  // split up event files - yes this is copy and pasted code
  ending = "neutron0_event.dat";
  start = eventfile.find(ending);
  if (start != string::npos)
    return eventfile.replace(start, ending.size(), "pulseid0.dat");

  ending = "neutron1_event.dat";
  start = eventfile.find(ending);
  if (start != string::npos)
    return eventfile.replace(start, ending.size(), "pulseid1.dat");

  return "";
}

203
//-----------------------------------------------------------------------------
204
static string generateMappingfileName(EventWorkspace_sptr &wksp) { //
205
  // get the name of the mapping file as set in the parameter files
206
207
  std::vector<string> temp =
      wksp->getInstrument()->getStringParameter("TS_mapping_file");
208
209
210
  if (temp.empty())
    return "";
  string mapping = temp[0];
211
212
213
214
  // Try to get it from the working directory
  Poco::File localmap(mapping);
  if (localmap.exists())
    return mapping;
215
216
217
218
219
220
221

  // Try to get it from the data directories
  string dataversion = Mantid::API::FileFinder::Instance().getFullPath(mapping);
  if (!dataversion.empty())
    return dataversion;

  // get a list of all proposal directories
222
  string instrument = wksp->getInstrument()->getName();
223
  Poco::File base("/SNS/" + instrument + "/");
224
  // try short instrument name
225
226
227
  if (!base.exists()) {
    instrument =
        Kernel::ConfigService::Instance().getInstrument(instrument).shortName();
Lynch, Vickie's avatar
Lynch, Vickie committed
228
    base = Poco::File("/SNS/" + instrument + "/");
229
230
231
    if (!base.exists())
      return "";
  }
232
233
234
235
236
  vector<string> dirs; // poco won't let me reuse temp
  base.list(dirs);

  // check all of the proposals for the mapping file in the canonical place
  const string CAL("_CAL");
237
  const size_t CAL_LEN = CAL.length(); // cache to make life easier
238
  vector<string> files;
239
240
241
242
  for (auto &dir : dirs) {
    if ((dir.length() > CAL_LEN) &&
        (dir.compare(dir.length() - CAL.length(), CAL.length(), CAL) == 0)) {
      if (Poco::File(base.path() + "/" + dir + "/calibrations/" + mapping)
243
              .exists())
244
        files.push_back(base.path() + "/" + dir + "/calibrations/" + mapping);
245
246
    }
  }
247
248
249

  if (files.empty())
    return "";
250
251
  else if (files.size() == 1)
    return files[0];
252
253
  else // just assume that the last one is the right one, this should never be
       // fired
254
    return *(files.rbegin());
255
256
}

257
258
259
260
namespace { // anonymous namespace
string getRunnumber(const string &filename) {
  // start by trimming the filename
  string runnumber(Poco::Path(filename).getBaseName());
261

262
263
264
265
  if (runnumber.find("neutron") >= string::npos)
    return "0";

  std::size_t left = runnumber.find("_");
266
  std::size_t right = runnumber.find("_", left + 1);
267

268
  return runnumber.substr(left + 1, right - left - 1);
269
270
}
}
271

272
//-----------------------------------------------------------------------------
273
/** Execute the algorithm */
274
void LoadEventPreNexus::exec() {
275
276
  // Check 'chunk' properties are valid, if set
  const int chunks = getProperty("TotalChunks");
277
  if (!isEmpty(chunks) && int(getProperty("ChunkNumber")) > chunks) {
278
279
280
    throw std::out_of_range("ChunkNumber cannot be larger than TotalChunks");
  }

281
  prog = new Progress(this, 0.0, 1.0, 100);
282

283
  // what spectra (pixel ID's) to load
284
  this->spectra_list = this->getProperty(PID_PARAM);
Nick Draper's avatar
re #100    
Nick Draper committed
285

286
287
  // the event file is needed in case the pulseid fileanme is empty
  string event_filename = this->getPropertyValue(EVENT_PARAM);
Nick Draper's avatar
re #100    
Nick Draper committed
288
  string pulseid_filename = this->getPropertyValue(PULSEID_PARAM);
289
  bool throwError = true;
290
  if (pulseid_filename.empty()) {
291
    pulseid_filename = generatePulseidName(event_filename);
292
293
294
295
    if (!pulseid_filename.empty()) {
      if (Poco::File(pulseid_filename).exists()) {
        this->g_log.information() << "Found pulseid file " << pulseid_filename
                                  << std::endl;
296
        throwError = false;
297
      } else {
298
        pulseid_filename = "";
299
      }
300
    }
301
302
  }

303
  prog->report("Loading Pulse ID file");
304
  this->readPulseidFile(pulseid_filename, throwError);
305

306
  this->openEventFile(event_filename);
Nick Draper's avatar
re #100    
Nick Draper committed
307

308
  prog->report("Creating output workspace");
Nick Draper's avatar
re #100    
Nick Draper committed
309
  // prep the output workspace
310
311
312
313
314
315
  EventWorkspace_sptr localWorkspace =
      EventWorkspace_sptr(new EventWorkspace());
  // Make sure to initialize.
  //   We can use dummy numbers for arguments, for event workspace it doesn't
  //   matter
  localWorkspace->initialize(1, 1, 1);
316

317
318
319
320
321
  // Set the units
  localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  localWorkspace->setYUnit("Counts");
  // TODO localWorkspace->setTitle(title);

322
323
  // Add the run_start property
  // Use the first pulse as the run_start time.
324
325
326
327
328
329
330
  if (this->num_pulses > 0) {
    // add the start of the run as a ISO8601 date/time string. The start = the
    // first pulse.
    // (this is used in LoadInstrument to find the right instrument file to
    // use).
    localWorkspace->mutableRun().addProperty(
        "run_start", pulsetimes[0].toISO8601String(), true);
331
332
  }

333
  // determine the run number and add it to the run object
334
335
  localWorkspace->mutableRun().addProperty("run_number",
                                           getRunnumber(event_filename));
336

337
  // Get the instrument!
338
  prog->report("Loading Instrument");
339
  this->runLoadInstrument(event_filename, localWorkspace);
340

341
  // load the mapping file
342
  prog->report("Loading Mapping File");
343
  string mapping_filename = this->getPropertyValue(MAP_PARAM);
344
345
346
  if (mapping_filename.empty()) {
    mapping_filename = generateMappingfileName(localWorkspace);
    if (!mapping_filename.empty())
347
348
      this->g_log.information() << "Found mapping file \"" << mapping_filename
                                << "\"" << std::endl;
349
  }
350
351
  this->loadPixelMap(mapping_filename);

352
  // Process the events into pixels
353
354
  this->procEvents(localWorkspace);

355
  // Save output
356
  this->setProperty<IEventWorkspace_sptr>(OUT_PARAM, localWorkspace);
357

358
359
  // Cleanup
  delete prog;
360
361
}

362
363
//-----------------------------------------------------------------------------
/** Load the instrument geometry File
364
 *  @param eventfilename :: Used to pick the instrument.
365
366
 *  @param localWorkspace :: MatrixWorkspace in which to put the instrument
 * geometry
367
 */
368
369
void LoadEventPreNexus::runLoadInstrument(const std::string &eventfilename,
                                          MatrixWorkspace_sptr localWorkspace) {
370
371
  // determine the instrument parameter file
  string instrument = Poco::Path(eventfilename).getFileName();
372
373
374
  size_t pos = instrument.rfind("_");   // get rid of 'event.dat'
  pos = instrument.rfind("_", pos - 1); // get rid of 'neutron'
  pos = instrument.rfind("_", pos - 1); // get rid of the run number
375
376
377
  instrument = instrument.substr(0, pos);

  // do the actual work
378
  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
379

380
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
381
  loadInst->setPropertyValue("InstrumentName", instrument);
382
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace);
383
  loadInst->setProperty("RewriteSpectraMap",
384
                        Mantid::Kernel::OptionalBool(false));
385
  loadInst->executeAsChildAlg();
386

387
388
  // Populate the instrument parameters in this workspace - this works around a
  // bug
389
  localWorkspace->populateInstrumentParameters();
390
391
}

392
//-----------------------------------------------------------------------------
393
394
395
/** Turn a pixel id into a "corrected" pixelid and period.
 *
 */
396
397
inline void LoadEventPreNexus::fixPixelId(PixelType &pixel,
                                          uint32_t &period) const {
398
  if (!this->using_mapping_file) { // nothing to do here
Nick Draper's avatar
re #100    
Nick Draper committed
399
400
401
402
403
404
405
406
407
    period = 0;
    return;
  }

  PixelType unmapped_pid = pixel % this->numpixel;
  period = (pixel - unmapped_pid) / this->numpixel;
  pixel = this->pixelmap[unmapped_pid];
}

408
//-----------------------------------------------------------------------------
409
/** Process the event file properly.
410
 * @param workspace :: EventWorkspace to write to.
411
 */
412
413
void LoadEventPreNexus::procEvents(
    DataObjects::EventWorkspace_sptr &workspace) {
Nick Draper's avatar
re #100    
Nick Draper committed
414
415
  this->num_error_events = 0;
  this->num_good_events = 0;
416
  this->num_ignored_events = 0;
Nick Draper's avatar
re #100    
Nick Draper committed
417

418
  // Default values in the case of no parallel
419
  size_t loadBlockSize = Mantid::Kernel::DEFAULT_BLOCK_SIZE * 2;
420
421
422

  shortest_tof = static_cast<double>(MAX_TOF_UINT32) * TOF_CONVERSION;
  longest_tof = 0.;
423

424
  // Initialize progress reporting.
425
  size_t numBlocks = (max_events + loadBlockSize - 1) / loadBlockSize;
426

427
  // We want to pad out empty pixels.
428
429
430
  detid2det_map detector_map;
  workspace->getInstrument()->getDetectors(detector_map);

431
432
433
434
435
436
  // -------------- Determine processing mode
  std::string procMode = getProperty("UseParallelProcessing");
  if (procMode == "Serial")
    parallelProcessing = false;
  else if (procMode == "Parallel")
    parallelProcessing = true;
437
438
439
440
441
442
443
  else {
    // Automatic determination. Loading serially (for me) is about 3 million
    // events per second,
    // (which is sped up by ~ x 3 with parallel processing, say 10 million per
    // second, e.g. 7 million events more per seconds).
    // compared to a setup time/merging time of about 10 seconds per million
    // detectors.
444
    double setUpTime = double(detector_map.size()) * 10e-6;
445
    parallelProcessing = ((double(max_events) / 7e6) > setUpTime);
446
447
    g_log.debug() << (parallelProcessing ? "Using" : "Not using")
                  << " parallel processing." << std::endl;
448
449
  }

450
451
  // determine maximum pixel id
  detid2det_map::iterator it;
452
  detid_max = 0; // seems like a safe lower bound
453
454
455
456
  for (it = detector_map.begin(); it != detector_map.end(); it++)
    if (it->first > detid_max)
      detid_max = it->first;

457
  // Pad all the pixels
458
  prog->report("Padding Pixels");
459
460
  this->pixel_to_wkspindex.reserve(
      detid_max + 1); // starting at zero up to and including detid_max
461
  // Set to zero
462
  this->pixel_to_wkspindex.assign(detid_max + 1, 0);
463
  size_t workspaceIndex = 0;
464
465
  for (it = detector_map.begin(); it != detector_map.end(); it++) {
    if (!it->second->isMonitor()) {
466
      this->pixel_to_wkspindex[it->first] = workspaceIndex;
467
      EventList &spec = workspace->getOrAddEventList(workspaceIndex);
468
      spec.addDetectorID(it->first);
469
      // Start the spectrum number at 1
470
      spec.setSpectrumNo(specid_t(workspaceIndex + 1));
471
      workspaceIndex += 1;
472
473
474
    }
  }

475
  // For slight speed up
476
477
  loadOnlySomeSpectra = (this->spectra_list.size() > 0);

478
  // Turn the spectra list into a map, for speed of access
Hahn, Steven's avatar
Hahn, Steven committed
479
480
  for (auto &spectrum : spectra_list)
    spectraLoadMap[spectrum] = true;
481

482
483
  CPUTimer tim;

484
485
  // --------------- Create the partial workspaces
  // ------------------------------------------
486
487
488
  // Vector of partial workspaces, for parallel processing.
  std::vector<EventWorkspace_sptr> partWorkspaces;
  std::vector<DasEvent *> buffers;
489
490

  /// Pointer to the vector of events
491
  typedef std::vector<TofEvent> *EventVector_pt;
492
  /// Bare array of arrays of pointers to the EventVectors
493
  EventVector_pt **eventVectors;
494
495
496

  /// How many threads will we use?
  size_t numThreads = 1;
497
  if (parallelProcessing)
498
499
500
501
502
503
    numThreads = size_t(PARALLEL_GET_MAX_THREADS);

  partWorkspaces.resize(numThreads);
  buffers.resize(numThreads);
  eventVectors = new EventVector_pt *[numThreads];

504
  // cppcheck-suppress syntaxError
505
  PRAGMA_OMP( parallel for if (parallelProcessing) )
506
  for (int i = 0; i < int(numThreads); i++) {
507
508
    // This is the partial workspace we are about to create (if in parallel)
    EventWorkspace_sptr partWS;
509
    if (parallelProcessing) {
510
511
      prog->report("Creating Partial Workspace");
      // Create a partial workspace
512
      partWS = EventWorkspace_sptr(new EventWorkspace());
513
514
515
516
      // Make sure to initialize.
      partWS->initialize(1, 1, 1);
      // Copy all the spectra numbers and stuff (no actual events to copy
      // though).
517
518
519
      partWS->copyDataFrom(*workspace);
      // Push it in the array
      partWorkspaces[i] = partWS;
520
    } else
521
522
      partWS = workspace;

523
    // Allocate the buffers
524
525
    buffers[i] = new DasEvent[loadBlockSize];

526
527
528
529
530
    // For each partial workspace, make an array where index = detector ID and
    // value = pointer to the events vector
    eventVectors[i] = new EventVector_pt[detid_max + 1];
    EventVector_pt *theseEventVectors = eventVectors[i];
    for (detid_t j = 0; j < detid_max + 1; j++) {
531
532
533
534
      size_t wi = pixel_to_wkspindex[j];
      // Save a POINTER to the vector<tofEvent>
      theseEventVectors[j] = &partWS->getEventList(wi).getEvents();
    }
535
536
  }

537
538
  g_log.debug() << tim << " to create " << partWorkspaces.size()
                << " workspaces for parallel loading." << std::endl;
539

540
  prog->resetNumSteps(numBlocks, 0.1, 0.8);
541
542
543

  // ---------------------------------- LOAD THE DATA --------------------------
  PRAGMA_OMP( parallel for schedule(dynamic, 1) if (parallelProcessing) )
544
  for (int blockNum = 0; blockNum < int(numBlocks); blockNum++) {
545
    PARALLEL_START_INTERUPT_REGION
546
547
548
549

    // Find the workspace for this particular thread
    EventWorkspace_sptr ws;
    size_t threadNum = 0;
550
    if (parallelProcessing) {
551
552
      threadNum = PARALLEL_THREAD_NUMBER;
      ws = partWorkspaces[threadNum];
553
    } else
554
      ws = workspace;
555

556
    // Get the buffer (for this thread)
557
    DasEvent *event_buffer = buffers[threadNum];
558

559
    // Get the speeding-up array of vector<tofEvent> where index = detid.
560
    EventVector_pt *theseEventVectors = eventVectors[threadNum];
561

562
    // Where to start in the file?
563
564
565
    size_t fileOffset = first_event + (loadBlockSize * blockNum);
    // May need to reduce size of last (or only) block
    size_t current_event_buffer_size =
566
567
568
        (blockNum == int(numBlocks - 1))
            ? (max_events - (numBlocks - 1) * loadBlockSize)
            : loadBlockSize;
569
570

    // Load this chunk of event data (critical block)
571
572
573
    PARALLEL_CRITICAL(LoadEventPreNexus_fileAccess) {
      current_event_buffer_size = eventfile->loadBlockAt(
          event_buffer, fileOffset, current_event_buffer_size);
574
575
576
    }

    // This processes the events. Can be done in parallel!
577
578
    procEventsLinear(ws, theseEventVectors, event_buffer,
                     current_event_buffer_size, fileOffset);
579
580
581

    // Report progress
    prog->report("Load Event PreNeXus");
582
583

    PARALLEL_END_INTERUPT_REGION
584
  }
585
  PARALLEL_CHECK_INTERUPT_REGION
586
587
  g_log.debug() << tim << " to load the data." << std::endl;

588
589
590
  // ---------------------------------- MERGE WORKSPACES BACK TOGETHER
  // --------------------------
  if (parallelProcessing) {
591
    PARALLEL_START_INTERUPT_REGION
592
    prog->resetNumSteps(workspace->getNumberHistograms(), 0.8, 0.95);
593
594
595

    size_t memoryCleared = 0;
    MemoryManager::Instance().releaseFreeMemory();
596

597
598
    // Merge all workspaces, index by index.
    PARALLEL_FOR_NO_WSP_CHECK()
599
    for (int iwi = 0; iwi < int(workspace->getNumberHistograms()); iwi++) {
600
601
602
      size_t wi = size_t(iwi);

      // The output event list.
603
      EventList &el = workspace->getEventList(wi);
604
605
606
607
      el.clear(false);

      // How many events will it have?
      size_t numEvents = 0;
608
      for (size_t i = 0; i < numThreads; i++)
609
610
611
        numEvents += partWorkspaces[i]->getEventList(wi).getNumberEvents();
      // This will avoid too much copying.
      el.reserve(numEvents);
612

613
      // Now merge the event lists
614
615
      for (size_t i = 0; i < numThreads; i++) {
        EventList &partEl = partWorkspaces[i]->getEventList(wi);
616
617
618
619
        el += partEl.getEvents();
        // Free up memory as you go along.
        partEl.clear(false);
      }
620

621
      // With TCMalloc, release memory when you accumulate enough to make sense
622
      PARALLEL_CRITICAL(LoadEventPreNexus_trackMemory) {
623
624
625
626
627
628
629
630
631
632
633
634
        memoryCleared += numEvents;
        if (memoryCleared > 10000000) // ten million events = about 160 MB
        {
          MemoryManager::Instance().releaseFreeMemory();
          memoryCleared = 0;
        }
      }
      prog->report("Merging Workspaces");
    }
    // Final memory release
    MemoryManager::Instance().releaseFreeMemory();
    g_log.debug() << tim << " to merge workspaces together." << std::endl;
635
    PARALLEL_END_INTERUPT_REGION
636
  }
637
  PARALLEL_CHECK_INTERUPT_REGION
638

639
  // Delete the buffers for each thread.
640
641
642
  for (size_t i = 0; i < numThreads; i++) {
    delete[] buffers[i];
    delete[] eventVectors[i];
643
  }
644
645
  delete[] eventVectors;
  // delete [] pulsetimes;
646

647
  prog->resetNumSteps(3, 0.94, 1.00);
648

649
  // finalize loading
650
  prog->report("Deleting Empty Lists");
651
  if (loadOnlySomeSpectra)
652
    workspace->deleteEmptyLists();
653

654
  prog->report("Setting proton charge");
655
  this->setProtonCharge(workspace);
656
  g_log.debug() << tim << " to set the proton charge log." << std::endl;
657

658
  // Make sure the MRU is cleared
659
  workspace->clearMRU();
660

661
  // Now, create a default X-vector for histogramming, with just 2 bins.
662
  Kernel::cow_ptr<MantidVec> axis;
663
  MantidVec &xRef = axis.access();
664
  xRef.resize(2);
665
  xRef[0] = shortest_tof - 1; // Just to make sure the bins hold it all
666
667
  xRef[1] = longest_tof + 1;
  workspace->setAllX(axis);
668
  this->pixel_to_wkspindex.clear();
669

670
  g_log.information() << "Read " << this->num_good_events << " events + "
671
672
673
674
                      << this->num_error_events << " errors"
                      << ". Shortest TOF: " << shortest_tof
                      << " microsec; longest TOF: " << longest_tof
                      << " microsec." << std::endl;
675
676
}

677
678
//-----------------------------------------------------------------------------
/** Linear-version of the procedure to process the event file properly.
679
 * @param workspace :: EventWorkspace to write to.
680
681
682
683
 * @param arrayOfVectors :: For speed up: this is an array, of size detid_max+1,
 * where the
 *        index is a pixel ID, and the value is a pointer to the
 * vector<tofEvent> in the given EventList.
684
685
686
 * @param event_buffer :: The buffer containing the DAS events
 * @param current_event_buffer_size :: The length of the given DAS buffer
 * @param fileOffset :: Value for an offset into the binary file
687
 */
688
689
690
691
692
void LoadEventPreNexus::procEventsLinear(
    DataObjects::EventWorkspace_sptr & /*workspace*/,
    std::vector<TofEvent> **arrayOfVectors, DasEvent *event_buffer,
    size_t current_event_buffer_size, size_t fileOffset) {
  // Starting pulse time
693
  DateAndTime pulsetime;
Peterson, Peter's avatar
Peterson, Peter committed
694
  int64_t pulse_i = 0;
695
  int64_t numPulses = static_cast<int64_t>(num_pulses);
696
697
698
  if (event_indices.size() < num_pulses) {
    g_log.warning()
        << "Event_indices vector is smaller than the pulsetimes array.\n";
Peterson, Peter's avatar
Peterson, Peter committed
699
    numPulses = static_cast<int64_t>(event_indices.size());
700
  }
701

702
703
704
  size_t local_num_error_events = 0;
  size_t local_num_ignored_events = 0;
  size_t local_num_good_events = 0;
705
706
  double local_shortest_tof =
      static_cast<double>(MAX_TOF_UINT32) * TOF_CONVERSION;
707
708
  double local_longest_tof = 0.;

709
  // process the individual events
710
711
  for (size_t i = 0; i < current_event_buffer_size; i++) {
    DasEvent &temp = *(event_buffer + i);
712
713
714
715
    PixelType pid = temp.pid;

    if ((pid & ERROR_PID) == ERROR_PID) // marked as bad
    {
716
      local_num_error_events++;
717
718
719
      continue;
    }

720
721
    // Covert the pixel ID from DAS pixel to our pixel ID
    if (this->using_mapping_file) {
722
723
724
      PixelType unmapped_pid = pid % this->numpixel;
      pid = this->pixelmap[unmapped_pid];
    }
725

726
    // Avoid segfaults for wrong pixel IDs
727
    if (pid > static_cast<PixelType>(detid_max)) {
728
729
730
731
      local_num_error_events++;
      continue;
    }

732
733
    // Now check if this pid we want to load.
    if (loadOnlySomeSpectra) {
Peterson, Peter's avatar
Peterson, Peter committed
734
      std::map<int64_t, bool>::iterator it;
735
      it = spectraLoadMap.find(pid);
736
737
      if (it == spectraLoadMap.end()) {
        // Pixel ID was not found, so the event is being ignored.
738
        local_num_ignored_events++;
739
740
741
742
743
        continue;
      }
    }

    // work with the good guys
744

745
746
747
    // Find the pulse time for this event index
    if (pulse_i < numPulses - 1) {
      // This is the total offset into the file
748
      size_t total_i = i + fileOffset;
749
750
751
752
      // Go through event_index until you find where the index increases to
      // encompass the current index. Your pulse = the one before.
      while (!((total_i >= event_indices[pulse_i]) &&
               (total_i < event_indices[pulse_i + 1]))) {
753
        pulse_i++;
754
        if (pulse_i >= (numPulses - 1))
755
756
          break;
      }
Janik Zikovsky's avatar
Janik Zikovsky committed
757

758
759
      // if (pulsetimes[pulse_i] != pulsetime)    std::cout << pulse_i << " at "
      // << pulsetimes[pulse_i] << "\n";
Janik Zikovsky's avatar
Janik Zikovsky committed
760

761
      // Save the pulse time at this index for creating those events
762
763
      pulsetime = pulsetimes[pulse_i];
    }
764

765
    double tof = static_cast<double>(temp.tof) * TOF_CONVERSION;
766
    TofEvent event(tof, pulsetime);
767

768
    // Find the overall max/min tof
769
770
771
772
    if (tof < local_shortest_tof)
      local_shortest_tof = tof;
    if (tof > local_longest_tof)
      local_longest_tof = tof;
773

774
775
776
    // The addEventQuickly method does not clear the cache, making things
    // slightly faster.
    // workspace->getEventList(this->pixel_to_wkspindex[pid]).addEventQuickly(event);
777

778
779
    // This is equivalent to
    // workspace->getEventList(this->pixel_to_wkspindex[pid]).addEventQuickly(event);
780
781
    // But should be faster as a bunch of these calls were cached.
    arrayOfVectors[pid]->push_back(event);
Janik Zikovsky's avatar
Janik Zikovsky committed
782

783
    // TODO work with period
784
    local_num_good_events++;
785

786
  } // for each event
787

788
  PARALLEL_CRITICAL(LoadEventPreNexus_global_statistics) {
789
790
791
792
793
794
795
796
    this->num_good_events += local_num_good_events;
    this->num_ignored_events += local_num_ignored_events;
    this->num_error_events += local_num_error_events;
    if (local_shortest_tof < shortest_tof)
      shortest_tof = local_shortest_tof;
    if (local_longest_tof > longest_tof)
      longest_tof = local_longest_tof;
  }
797
798
799
800
}

//-----------------------------------------------------------------------------
/// Comparator for sorting dasevent lists
801
bool intermediatePixelIDComp(IntermediateEvent x, IntermediateEvent y) {
802
803
804
  return (x.pid < y.pid);
}

805
//-----------------------------------------------------------------------------
806
/**
807
808
809
810
 * Add a sample environment log for the proton chage (charge of the pulse in
 *picoCoulombs)
 * and set the scalar value (total proton charge, microAmps*hours, on the
 *sample)
811
 *
812
 * @param workspace :: Event workspace to set the proton charge on
813
 */
814
815
void LoadEventPreNexus::setProtonCharge(
    DataObjects::EventWorkspace_sptr &workspace) {
816
  if (this->proton_charge.empty()) // nothing to do
817
    return;
818

819
  Run &run = workspace->mutableRun();
820

821
822
823
  // Add the proton charge entries.
  TimeSeriesProperty<double> *log =
      new TimeSeriesProperty<double>("proton_charge");
824
  log->setUnits("picoCoulombs");
825

826