FilterEventsByLogValuePreNexus.cpp 89.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include "MantidDataHandling/FilterEventsByLogValuePreNexus.h"
#include <algorithm>
#include <sstream>
#include <stdexcept>
#include <functional>
#include <iostream>
#include <set>
#include <vector>
#include <Poco/File.h>
#include <Poco/Path.h>
#include <boost/timer.hpp>
#include "MantidAPI/FileFinder.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/FileValidator.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/Glob.h"
#include "MantidAPI/FileProperty.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/BinaryFile.h"
#include "MantidKernel/System.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidGeometry/IDetector.h"
#include "MantidKernel/CPUTimer.h"
#include "MantidKernel/VisibleWhenProperty.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"

#include <algorithm>
#include <sstream>
#include "MantidAPI/MemoryManager.h"

Zhou, Wenduo's avatar
Zhou, Wenduo committed
41
42
// #define DBOUT

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
namespace Mantid {
namespace DataHandling {

DECLARE_FILELOADER_ALGORITHM(FilterEventsByLogValuePreNexus)

using namespace Kernel;
using namespace API;
using namespace DataObjects;
using namespace Geometry;
using boost::posix_time::ptime;
using boost::posix_time::time_duration;
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;
using std::string;
using std::vector;

//----------------------------------------------------------------------------------------------
// constants for locating the parameters to use in execution
//----------------------------------------------------------------------------------------------
static const string EVENT_PARAM("EventFilename");
static const string PULSEID_PARAM("PulseidFilename");
static const string MAP_PARAM("MappingFilename");
static const string PID_PARAM("SpectrumList");
static const string PARALLEL_PARAM("UseParallelProcessing");
static const string BLOCK_SIZE_PARAM("LoadingBlockSize");
static const string OUT_PARAM("OutputWorkspace");
/// All pixel ids with matching this mask are errors.
static const PixelType ERROR_PID = 0x80000000;
/// The maximum possible tof as native type
static const uint32_t MAX_TOF_UINT32 = std::numeric_limits<uint32_t>::max();
/// 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.;
/// Veto mask as 0xFF000000000
static const uint64_t VETOFLAG(72057594037927935);

static const string EVENT_EXTS[] = {
    "_neutron_event.dat",  "_neutron0_event.dat", "_neutron1_event.dat",
    "_neutron2_event.dat", "_neutron3_event.dat", "_live_neutron_event.dat"};
static const string PULSE_EXTS[] = {"_pulseid.dat",  "_pulseid0.dat",
                                    "_pulseid1.dat", "_pulseid2.dat",
                                    "_pulseid3.dat", "_live_pulseid.dat"};
static const int NUM_EXT = 6;

//----------------------------------------------------------------------------------------------
// Functions to deal with file name and run information
//----------------------------------------------------------------------------------------------

//----------------------------------------------------------------------------------------------
/** Get run number
  */
static string getRunnumber(const string &filename) {
  // start by trimming the filename
  string runnumber(Poco::Path(filename).getBaseName());
105

106
107
  if (runnumber.find("neutron") >= string::npos)
    return "0";
108

109
110
111
112
113
  std::size_t left = runnumber.find("_");
  std::size_t right = runnumber.find("_", left + 1);

  return runnumber.substr(left + 1, right - left - 1);
}
114

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
//----------------------------------------------------------------------------------------------
/** Generate pulse ID
  */
static string generatePulseidName(string eventfile) {
  // initialize vector of endings and put live at the beginning
  vector<string> eventExts(EVENT_EXTS, EVENT_EXTS + NUM_EXT);
  std::reverse(eventExts.begin(), eventExts.end());
  vector<string> pulseExts(PULSE_EXTS, PULSE_EXTS + NUM_EXT);
  std::reverse(pulseExts.begin(), pulseExts.end());

  // look for the correct ending
  for (std::size_t i = 0; i < eventExts.size(); ++i) {
    size_t start = eventfile.find(eventExts[i]);
    if (start != string::npos)
      return eventfile.replace(start, eventExts[i].size(), pulseExts[i]);
130
131
  }

132
133
134
  // give up and return nothing
  return "";
}
135

136
137
138
139
140
141
142
143
//----------------------------------------------------------------------------------------------
/** Generate mapping file name
  */
static string generateMappingfileName(EventWorkspace_sptr &wksp) {
  // get the name of the mapping file as set in the parameter files
  std::vector<string> temp =
      wksp->getInstrument()->getStringParameter("TS_mapping_file");
  if (temp.empty())
144
145
    return "";

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
  string mapping = temp[0];
  // Try to get it from the working directory
  Poco::File localmap(mapping);
  if (localmap.exists())
    return mapping;

  // 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
  string instrument = wksp->getInstrument()->getName();
  Poco::File base("/SNS/" + instrument + "/");
  // try short instrument name
  if (!base.exists()) {
    return "";
163
#if 0
164
165
166
167
      instrument = Kernel::ConfigService::Instance().getInstrument(instrument).shortName();
      base = Poco::File("/SNS/" + instrument + "/");
      if (!base.exists())
        return "";
168
#endif
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  }
  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");
  const size_t CAL_LEN = CAL.length(); // cache to make life easier
  vector<string> files;
  for (size_t i = 0; i < dirs.size(); ++i) {
    if ((dirs[i].length() > CAL_LEN) &&
        (dirs[i].compare(dirs[i].length() - CAL.length(), CAL.length(), CAL) ==
         0)) {
      if (Poco::File(base.path() + "/" + dirs[i] + "/calibrations/" + mapping)
              .exists())
        files.push_back(base.path() + "/" + dirs[i] + "/calibrations/" +
                        mapping);
185
186
187
    }
  }

188
189
190
191
192
193
194
195
  if (files.empty())
    return "";
  else if (files.size() == 1)
    return files[0];
  else // just assume that the last one is the right one, this should never be
       // fired
    return *(files.rbegin());
}
196

197
198
199
//----------------------------------------------------------------------------------------------
// Member functions
//----------------------------------------------------------------------------------------------
200

201
202
203
204
//----------------------------------------------------------------------------------------------
/** Constructor
*/
FilterEventsByLogValuePreNexus::FilterEventsByLogValuePreNexus()
205
    : Mantid::API::IFileLoader<Kernel::FileDescriptor>(), proton_charge_tot(0),
206
      detid_max(0), eventfile(NULL), num_events(0), num_pulses(0), numpixel(0),
207
208
209
210
211
212
213
      num_good_events(0), num_error_events(0), num_bad_events(0),
      num_wrongdetid_events(0), num_ignored_events(0), first_event(0),
      m_maxNumEvents(0), using_mapping_file(false), loadOnlySomeSpectra(false),
      longest_tof(0.0), shortest_tof(0.0), parallelProcessing(false),
      pulsetimesincreasing(false), m_throwError(true), m_examEventLog(false),
      m_pixelid2exam(0), m_numevents2write(0), m_freqHz(0), istep(0),
      m_dbPixelID(0), m_useDBOutput(false), m_corretctTOF(false) {}
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

//----------------------------------------------------------------------------------------------
/** Desctructor
 */
FilterEventsByLogValuePreNexus::~FilterEventsByLogValuePreNexus() {
  delete this->eventfile;
}

//----------------------------------------------------------------------------------------------
/** Return the confidence with with this algorithm can load the file
 *  @param descriptor A descriptor for the file
 *  @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int FilterEventsByLogValuePreNexus::confidence(
    Kernel::FileDescriptor &descriptor) const {
  if (descriptor.extension().rfind("dat") == std::string::npos)
    return 0;

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

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

  if (filesize % objSize == 0)
    return 10;
  else
    return 0;
}

//----------------------------------------------------------------------------------------------
/**  Initialize the algorithm
254
  */
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
void FilterEventsByLogValuePreNexus::init() {
  // File files to use
  vector<string> eventExts(EVENT_EXTS, EVENT_EXTS + NUM_EXT);
  declareProperty(
      new FileProperty(EVENT_PARAM, "", FileProperty::Load, eventExts),
      "The name of the neutron event file to read, including its full or "
      "relative path. In most cases, the file typically ends in "
      "neutron_event.dat (N.B. case sensitive if running on Linux).");
  vector<string> pulseExts(PULSE_EXTS, PULSE_EXTS + NUM_EXT);
  declareProperty(new FileProperty(PULSEID_PARAM, "",
                                   FileProperty::OptionalLoad, pulseExts),
                  "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.");

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

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

  // Loading option
  std::vector<std::string> propOptions;
  propOptions.push_back("Auto");
  propOptions.push_back("Serial");
  propOptions.push_back("Parallel");
  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.");

  // the output workspace name
  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]].");

  // Optional output table workspace
  declareProperty(new WorkspaceProperty<ITableWorkspace>(
                      "EventLogTableWorkspace", "", PropertyMode::Optional),
                  "Optional output table workspace containing the event log "
                  "(pixel) information. ");

  //
  std::vector<std::string> vecfunmode;
  vecfunmode.push_back("LoadData");
  vecfunmode.push_back("Filter");
  vecfunmode.push_back("ExamineEventLog");
  declareProperty("FunctionMode", "LoadData",
                  boost::make_shared<StringListValidator>(vecfunmode),
                  "Function mode for different purpose. ");

  declareProperty("PixelIDtoExamine", EMPTY_INT(),
                  "Pixel ID for the events to be examined. ");

  declareProperty("NumberOfEventsToExamine", EMPTY_INT(),
                  "Number of events on the pixel ID to get examined. ");

  declareProperty(new ArrayProperty<int>("LogPixelIDs"),
                  "Pixel IDs for event log. Must have 2 (or more) entries. ");

  declareProperty(
      new ArrayProperty<std::string>("LogPIxelTags"),
      "Pixel ID tags for event log. Must have same items as 'LogPixelIDs'. ");

  declareProperty("AcceleratorFrequency", 60, "Freuqency of the accelerator at "
                                              "which the experiment runs. It "
                                              "can 20, 30 or 60.");

  declareProperty("CorrectTOFtoSample", false,
                  "Correct TOF to sample position. ");

  declareProperty("DBPixelID", EMPTY_INT(),
                  "ID of the pixel (detector) for debug output. ");

  return;
}
352

353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
//----------------------------------------------------------------------------------------------
/** Execute the algorithm
 * 1. check all the inputs
 * 2. create an EventWorkspace object
 * 3. process events
 * 4. set out output
 */
void FilterEventsByLogValuePreNexus::exec() {
  // Process inputs
  prog = new Progress(this, 0.0, 1.0, 100);
  processProperties();

  // Read input files
  prog->report("Loading Pulse ID file");
  readPulseidFile(m_pulseIDFileName, m_throwError);

  prog->report("Loading Event File");
  openEventFile(m_eventFileName);

  // Correct wrong event index in loaded eventindexes
  unmaskVetoEventIndexes();

  // Find out the frequency of the frequency
  int runfreq = findRunFrequency();
  if (m_freqHz != runfreq) {
    if (m_freqHz % runfreq == 0) {
      int frame = m_freqHz / runfreq;
      g_log.warning() << "Input frequency " << m_freqHz
                      << " is different from data. "
                      << "It is forced to use input frequency, while all "
                         "events' pulse time will be "
                      << "set to " << frame << "-th freme. "
                      << "\n";
    } else {
      throw std::runtime_error("Operation frequency is not self-consistent");
    }
389
  }
390
391
392
393
394
395
396
397
398
399
400
401
  istep = 60 / m_freqHz;

  // Create and set up output EventWorkspace
  localWorkspace = setupOutputEventWorkspace();
  if (m_functionMode == "Filter")
    localWorkspaceBA = setupOutputEventWorkspace();

  // Process the events into pixels
  if (m_functionMode == "Filter") {
    filterEvents();
  } else {
    procEvents(localWorkspace);
402
403
  }

404
405
406
407
408
409
410
411
412
413
414
415
  // set that the sort order on the event lists
  if (this->num_pulses > 0 && this->pulsetimesincreasing) {
    const int64_t numberOfSpectra = localWorkspace->getNumberHistograms();
    PARALLEL_FOR_NO_WSP_CHECK()
    for (int64_t i = 0; i < numberOfSpectra; i++) {
      PARALLEL_START_INTERUPT_REGION
      localWorkspace->getEventListPtr(i)
          ->setSortOrder(DataObjects::PULSETIME_SORT);
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }
416

417
418
419
420
421
422
423
424
  // Save output
  setProperty<IEventWorkspace_sptr>(OUT_PARAM, localWorkspace);
  if (m_functionMode == "Filter") {
    declareProperty(new WorkspaceProperty<IEventWorkspace>(
                        "OutputFilteredWorkspace", "WS_A", Direction::Output),
                    "");
    setProperty<IEventWorkspace_sptr>("OutputFilteredWorkspace",
                                      localWorkspaceBA);
425
426
  }

427
428
  // Add fast frequency sample environment (events) data to workspace's log
  processEventLogs();
429

430
431
  // -1. Cleanup
  delete prog;
432

433
434
  return;
} // exec()
435

436
437
438
439
440
441
442
443
444
445
//----------------------------------------------------------------------------------------------
/** Process input properties
  */
void FilterEventsByLogValuePreNexus::processProperties() {
  // Process and check input properties
  // Check 'chunk' properties are valid, if set
  const int chunks = getProperty("TotalChunks");
  if (!isEmpty(chunks) && int(getProperty("ChunkNumber")) > chunks) {
    throw std::out_of_range("ChunkNumber cannot be larger than TotalChunks");
  }
446

447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
  // What spectra (pixel ID's) to load
  this->spectra_list = this->getProperty(PID_PARAM);

  // The event file is needed in case the pulseid fileanme is empty
  m_eventFileName = this->getPropertyValue(EVENT_PARAM);

  // Pulse ID file
  m_pulseIDFileName = this->getPropertyValue(PULSEID_PARAM);
  m_throwError = true;

  if (m_pulseIDFileName.empty()) {
    // Pulse ID file is not given: generate by routine
    m_pulseIDFileName = generatePulseidName(m_eventFileName);
    if (!m_pulseIDFileName.empty()) {
      // Check existence of pulse ID file with generated name
      if (Poco::File(m_pulseIDFileName).exists()) {
        g_log.information() << "Found pulseid file " << m_pulseIDFileName
                            << "\n";
        m_throwError = false;
      } else {
        m_pulseIDFileName = "";
        g_log.warning("Generated pulse ID file name does not point to an "
                      "existing file. ");
470
      }
471
472
    } else {
      g_log.warning("Generated an empty pulse ID file. ");
473
    }
474
  }
475

476
  m_functionMode = getPropertyValue("FunctionMode");
477

478
479
  m_pixelid2exam = getProperty("PixelIDtoExamine");
  m_numevents2write = getProperty("NumberOfEventsToExamine");
480

481
482
483
484
485
486
  // Check whether option function mode is valid
  m_examEventLog = false;
  if (m_functionMode == "ExamineEventLog") {
    bool nogo = false;
    if (isEmpty(m_pixelid2exam)) {
      nogo = true;
487
488
    }

489
490
491
492
493
    if (nogo) {
      g_log.warning()
          << "In functional mode ExamineEventLog, pixel ID must be given!"
          << "\n";
      throw std::runtime_error("Incorrect input.");
494
495
    }

496
497
498
499
500
501
502
503
504
505
506
507
    m_examEventLog = true;
  } else if (m_functionMode == "Filter") {
    m_vecLogPixelID = getProperty("LogPixelIDs");
    m_vecLogPixelTag = getProperty("LogPIxelTags");

    if (m_vecLogPixelID.size() < 2) {
      throw std::runtime_error(
          "Input log pixel IDs must have more than 2 entries. ");
    } else if (m_vecLogPixelID.size() != m_vecLogPixelTag.size()) {
      throw std::runtime_error(
          "Input log pixel tags must have the same number of items as "
          "log pixe IDs. ");
508
    }
509
  }
510

511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
  //---------------------------------------------------------------------------
  // Load partial spectra
  //---------------------------------------------------------------------------
  // For slight speed up
  loadOnlySomeSpectra = (this->spectra_list.size() > 0);

  // Turn the spectra list into a map, for speed of access
  for (std::vector<int64_t>::iterator it = spectra_list.begin();
       it != spectra_list.end(); it++)
    spectraLoadMap[*it] = true;

  //---------------------------------------------------------------------------
  // Other features
  //---------------------------------------------------------------------------
  // Accelerator frquency
  m_freqHz = getProperty("AcceleratorFrequency");
  if (m_freqHz != 20 && m_freqHz != 30 && m_freqHz != 60)
    throw runtime_error("Only 20, 30 and 60Hz are supported. ");

  int tempint = getProperty("DBPixelID");
  if (isEmpty(tempint))
    m_useDBOutput = false;
  else {
    m_useDBOutput = true;
    m_dbPixelID = static_cast<int64_t>(tempint);
  }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
537

538
  m_corretctTOF = getProperty("CorrectTOFtoSample");
Zhou, Wenduo's avatar
Zhou, Wenduo committed
539

540
541
  return;
} // END of processProperties
542

543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
//----------------------------------------------------------------------------------------------
/** Create, initialize and set up output EventWorkspace
  */
DataObjects::EventWorkspace_sptr
FilterEventsByLogValuePreNexus::setupOutputEventWorkspace() {
  // Create and initialize output EventWorkspace
  prog->report("Creating output workspace");

  EventWorkspace_sptr tempworkspace;

  tempworkspace = EventWorkspace_sptr(new EventWorkspace());
  // Make sure to initialize. We can use dummy numbers for arguments, for event
  // workspace it doesn't matter
  tempworkspace->initialize(1, 1, 1);
  // Set the units and title
  tempworkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  tempworkspace->setYUnit("Counts");
  tempworkspace->setTitle("Dummy Title");

  // Add some properties to output workspace, including
  //   the run_start property (Use the first pulse as the run_start time)
  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).
Zhou, Wenduo's avatar
Zhou, Wenduo committed
568
#if 0
569
570
571
      DateAndTime pulse0 = pulsetimes[0];
      g_log.notice() << "Pulse time 0 = " <<  pulse0.totalNanoseconds() << "\n";
#endif
572
573
574
    tempworkspace->mutableRun().addProperty(
        "run_start", pulsetimes[0].toISO8601String(), true);
  }
575

576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
  //   the run number and add it to the run object
  tempworkspace->mutableRun().addProperty("run_number",
                                          getRunnumber(m_eventFileName));

  // Add the instrument!
  prog->report("Loading Instrument");
  this->runLoadInstrument(m_eventFileName, tempworkspace);

  // Load the mapping file
  prog->report("Loading Mapping File");
  string mapping_filename = this->getPropertyValue(MAP_PARAM);
  if (mapping_filename.empty()) {
    // No mapping file given: genrate mapping file name by routine
    mapping_filename = generateMappingfileName(tempworkspace);
    if (!mapping_filename.empty())
      g_log.information() << "Found mapping file \"" << mapping_filename << "\""
                          << "\n";
    else
      g_log.warning("No mapping file is generated. ");
595
  }
596
597
  // if (!mapping_filename.empty())
  loadPixelMap(mapping_filename);
598

599
600
  return tempworkspace;
}
601

602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
//----------------------------------------------------------------------------------------------
/** Process imbed logs (marked by bad pixel IDs)
  * (1) Add special event log to workspace log
  * (2) (Optionally) do statistic to each pixel
  * (3) (Optionally) write out information
  */
void FilterEventsByLogValuePreNexus::processEventLogs() {
  std::set<PixelType>::iterator pit;
  std::map<PixelType, size_t>::iterator mit;
  for (pit = this->wrongdetids.begin(); pit != this->wrongdetids.end(); ++pit) {
    // Convert Pixel ID to 'wrong detectors ID' map's index
    PixelType pid = *pit;
    mit = this->wrongdetidmap.find(pid);
    size_t mindex = mit->second;
    if (mindex > this->wrongdetid_pulsetimes.size()) {
      g_log.error() << "Wrong Index " << mindex << " for Pixel " << pid
                    << std::endl;
      throw std::invalid_argument("Wrong array index for pixel from map");
    } else {
      g_log.information() << "Processing imbed log marked by Pixel " << pid
                          << " with size = "
                          << this->wrongdetid_pulsetimes[mindex].size()
                          << std::endl;
    }
626

627
628
629
630
    // Generate the log name
    std::stringstream ssname;
    ssname << "Pixel" << pid;
    std::string logname = ssname.str();
631

632
633
    // Add this map entry to log
    addToWorkspaceLog(logname, mindex);
634

635
636
    // Do some statistic to this event log
    doStatToEventLog(mindex);
637

638
639
    g_log.information() << "Added Log " << logname
                        << " to output workspace. \n";
640

641
  } // ENDFOR pit
642

643
644
645
646
647
648
649
  // Output table workspace
  std::string evlog = getPropertyValue("EventLogTableWorkspace");
  if (!evlog.empty()) {
    // Initialize table workspace
    TableWorkspace_sptr evtablews = boost::make_shared<TableWorkspace>();
    evtablews->addColumn("int", "Pixel-ID");
    evtablews->addColumn("int", "NumberOfEvents");
650

651
652
653
654
655
656
    // Add information rows
    std::map<PixelType, size_t>::iterator git;
    for (git = this->wrongdetidmap.begin(); git != this->wrongdetidmap.end();
         ++git) {
      PixelType tmpid = git->first;
      size_t vindex = git->second;
657

658
659
660
      TableRow temprow = evtablews->appendRow();
      temprow << static_cast<int>(tmpid)
              << static_cast<int>(wrongdetid_pulsetimes[vindex].size());
661
662
    }

663
664
665
    // Set property
    setProperty("EventLogTableWorkspace",
                boost::dynamic_pointer_cast<ITableWorkspace>(evtablews));
666
667
  }

668
669
  return;
}
670

671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
//----------------------------------------------------------------------------------------------
/** Add absolute time series to log
  * @param logtitle :: title of the log to be inserted to workspace
  * @param mindex ::  index of the the series in the wrong detectors map
  */
void FilterEventsByLogValuePreNexus::addToWorkspaceLog(std::string logtitle,
                                                       size_t mindex) {
  // Create TimeSeriesProperty
  TimeSeriesProperty<double> *property =
      new TimeSeriesProperty<double>(logtitle);

  // Add entries
  size_t nbins = this->wrongdetid_pulsetimes[mindex].size();
  for (size_t k = 0; k < nbins; k++) {
    double tof = this->wrongdetid_tofs[mindex][k];
    DateAndTime pulsetime = wrongdetid_pulsetimes[mindex][k];
    int64_t abstime_ns =
        pulsetime.totalNanoseconds() + static_cast<int64_t>(tof * 1000);
    DateAndTime abstime(abstime_ns);

    double value = tof;

    property->addValue(abstime, value);
  } // ENDFOR

  // Add property to workspace
  localWorkspace->mutableRun().addProperty(property, false);

  g_log.information() << "Size of Property " << property->name() << " = "
                      << property->size() << " vs Original Log Size = " << nbins
                      << "\n";
702

703
704
  return;
}
705

706
707
708
709
710
711
712
713
714
715
716
//----------------------------------------------------------------------------------------------
/** Perform statistics to event (wrong pixel ID) logs
  * @param mindex ::  index of the the series in the list
  */
void FilterEventsByLogValuePreNexus::doStatToEventLog(size_t mindex) {
  // Create a vector of event log time entries
  size_t nbins = this->wrongdetid_pulsetimes[mindex].size();
  if (nbins <= 2) {
    g_log.warning() << "Event log of map index " << mindex << " has " << nbins
                    << " entries.  There is no need to do statistic on it. "
                    << "\n";
717
718
  }

719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
  std::vector<int64_t> vec_logtime(nbins);
  for (size_t i = 0; i < nbins; ++i) {
    DateAndTime ptime = wrongdetid_pulsetimes[mindex][i];
    int64_t templogtime =
        ptime.totalNanoseconds() +
        static_cast<int64_t>(wrongdetid_tofs[mindex][i] * 1000.);

    vec_logtime[i] = templogtime;
  } // ENDFOR

  // Sort
  std::sort(vec_logtime.begin(), vec_logtime.end());

  // Do statistic
  int64_t min_dt = vec_logtime[1] - vec_logtime[0];
  int64_t max_dt = min_dt;
  int64_t sum_dt = min_dt;
  int64_t numzeros = 0;
  for (size_t i = 2; i < nbins; ++i) {
    int64_t temp_dt = vec_logtime[i] - vec_logtime[i - 1];
    if (temp_dt == 0)
      ++numzeros;

    sum_dt += temp_dt;

    if (temp_dt < min_dt)
      min_dt = temp_dt;
    else if (temp_dt > max_dt)
      max_dt = temp_dt;
748
749
  }

750
751
752
753
754
755
756
757
758
  if ( nbins - 1 ) {
    double avg_dt = static_cast<double>(sum_dt) / static_cast<double>(nbins - 1);
    g_log.information() << "Event log of map index " << mindex
      << ": Avg(dt) = " << avg_dt * 1.0E-9
      << ", Min(dt) = " << static_cast<double>(min_dt) * 1.0E-9
      << ", Max(dt) = " << static_cast<double>(max_dt) * 1.0E-9
      << "\n";
  } else {
    g_log.information() << "Event log of map index " << mindex
759
      << ": Avg(dt) = " << static_cast<double>(sum_dt) * 1.0E-9
760
761
762
763
764
      << ", Min(dt) = " << static_cast<double>(min_dt) * 1.0E-9
      << ", Max(dt) = " << static_cast<double>(max_dt) * 1.0E-9
      << "\n";

  }
765

766
767
  g_log.information() << "Number of zero-interval eveng log = " << numzeros
                      << "\n";
768

769
770
  return;
}
771

772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
//-----------------------------------------------------------------------------
/** Load the instrument geometry File
 *  @param eventfilename :: Used to pick the instrument.
 *  @param localWorkspace :: MatrixWorkspace in which to put the instrument
 * geometry
 */
void FilterEventsByLogValuePreNexus::runLoadInstrument(
    const std::string &eventfilename, MatrixWorkspace_sptr localWorkspace) {
  // start by getting just the filename
  string instrument = Poco::Path(eventfilename).getFileName();

  // initialize vector of endings and put live at the beginning
  vector<string> eventExts(EVENT_EXTS, EVENT_EXTS + NUM_EXT);
  std::reverse(eventExts.begin(), eventExts.end());

  for (size_t i = 0; i < eventExts.size(); ++i) {
    size_t pos = instrument.find(eventExts[i]);
    if (pos != string::npos) {
      instrument = instrument.substr(0, pos);
      break;
    }
  }
794

795
796
797
  // determine the instrument parameter file
  size_t pos = instrument.rfind("_"); // get rid of the run number
  instrument = instrument.substr(0, pos);
798

799
800
  // do the actual work
  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
801

802
803
804
805
806
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  loadInst->setPropertyValue("InstrumentName", instrument);
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace);
  loadInst->setProperty("RewriteSpectraMap", false);
  loadInst->executeAsChildAlg();
807

808
809
810
811
  // Populate the instrument parameters in this workspace - this works around a
  // bug
  localWorkspace->populateInstrumentParameters();
}
812

813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
//----------------------------------------------------------------------------------------------
/** Process the event file properly.
 * @param workspace :: EventWorkspace to write to.
 */
void FilterEventsByLogValuePreNexus::procEvents(
    DataObjects::EventWorkspace_sptr &workspace) {
  // Initialize stat parameters
  this->num_error_events = 0;
  this->num_good_events = 0;
  this->num_ignored_events = 0;
  this->num_bad_events = 0;
  this->num_wrongdetid_events = 0;

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

  // -------------------------------------------------------------------
  // Set up instrument related parameters such as detector map and etc.
  // We want to pad out empty pixels.
  //--------------------------------------------------------------------
  detid2det_map detector_map;
  workspace->getInstrument()->getDetectors(detector_map);

  // Determine maximum pixel id
  detid2det_map::iterator it;
  detid_max = 0; // seems like a safe lower bound
  for (it = detector_map.begin(); it != detector_map.end(); it++)
    if (it->first > detid_max)
      detid_max = it->first;

  // Pad all the pixels
  prog->report("Padding Pixels");
  this->pixel_to_wkspindex.reserve(
      detid_max + 1); // starting at zero up to and including detid_max
  // Set to zero
  this->pixel_to_wkspindex.assign(detid_max + 1, 0);
  size_t workspaceIndex = 0;
  for (it = detector_map.begin(); it != detector_map.end(); it++) {
    if (!it->second->isMonitor()) {
      // Add non-monitor detector ID
      this->pixel_to_wkspindex[it->first] = workspaceIndex;
      EventList &spec = workspace->getOrAddEventList(workspaceIndex);
      spec.addDetectorID(it->first);
      // Start the spectrum number at 1
      spec.setSpectrumNo(specid_t(workspaceIndex + 1));
      workspaceIndex += 1;
859
    }
860
  }
861

862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
  // ----------------------------------------------------------------
  // Determine processing mode and file-loading parameters
  //------------------------------------------------------------------
  // Set up some default values in the case of no parallel
  size_t loadBlockSize = Mantid::Kernel::DEFAULT_BLOCK_SIZE * 2;
  size_t numBlocks = (m_maxNumEvents + loadBlockSize - 1) / loadBlockSize;

  std::string procMode = getProperty("UseParallelProcessing");
  if (procMode == "Serial") {
    parallelProcessing = false;
  } else if (procMode == "Parallel") {
    parallelProcessing = true;
  } 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.
    double setUpTime = double(detector_map.size()) * 10e-6;
    parallelProcessing = ((double(m_maxNumEvents) / 7e6) > setUpTime);
    g_log.information() << (parallelProcessing ? "Using" : "Not using")
                        << " parallel processing."
                        << "\n";
  }
887

888
889
890
891
892
  if (m_functionMode == "ExamineEventLog" && parallelProcessing) {
    parallelProcessing = false;
    g_log.notice("In function mode 'ExamineEventLog', processing mode is "
                 "forced to serial. ");
  }
893
894
895
896
897
898
899
900
901
902

#if 0
    //For slight speed up
    loadOnlySomeSpectra = (this->spectra_list.size() > 0);

    //Turn the spectra list into a map, for speed of access
    for (std::vector<int64_t>::iterator it = spectra_list.begin(); it != spectra_list.end(); it++)
      spectraLoadMap[*it] = true;
#endif

903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
  CPUTimer tim;

  // -------------------------------------------------------------------
  // Create the partial workspaces
  //--------------------------------------------------------------------
  // Vector of partial workspaces, for parallel processing.
  std::vector<EventWorkspace_sptr> partWorkspaces;
  std::vector<DasEvent *> buffers;

  /// Pointer to the vector of events
  typedef std::vector<TofEvent> *EventVector_pt;
  /// Bare array of arrays of pointers to the EventVectors
  EventVector_pt **eventVectors;

  /// How many threads will we use?
  size_t numThreads = 1;
  if (parallelProcessing)
    numThreads = size_t(PARALLEL_GET_MAX_THREADS);

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

  // Processing by number of threads
  g_log.information() << "Processing input event preNexus by " << numThreads
                      << " threads"
                      << " in " << numBlocks << " blocks. "
                      << "\n";
931

932
  // cppcheck-suppress syntaxError
933
    PRAGMA_OMP( parallel for schedule(dynamic, 1) if (parallelProcessing) )
934
    for (int i = 0; i < int(numThreads); i++) {
935
936
937
      // This is the partial workspace we are about to create (if in parallel)
      EventWorkspace_sptr partWS;

938
      if (parallelProcessing) {
939
940
941
        prog->report("Creating Partial Workspace");
        // Create a partial workspace
        partWS = EventWorkspace_sptr(new EventWorkspace());
942
943
944
945
        // Make sure to initialize.
        partWS->initialize(1, 1, 1);
        // Copy all the spectra numbers and stuff (no actual events to copy
        // though).
946
947
948
        partWS->copyDataFrom(*workspace);
        // Push it in the array
        partWorkspaces[i] = partWS;
949
      } else
950
951
        partWS = workspace;

952
      // Allocate the buffers
953
954
      buffers[i] = new DasEvent[loadBlockSize];

955
956
957
958
959
      // 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++) {
960
961
962
963
964
965
        size_t wi = pixel_to_wkspindex[j];
        // Save a POINTER to the vector<tofEvent>
        theseEventVectors[j] = &partWS->getEventList(wi).getEvents();
      }
    } // END FOR [Threads]

966
967
    g_log.information() << tim << " to create " << partWorkspaces.size()
                        << " workspaces for parallel loading."
968
969
                        << "\n";

970
    prog->resetNumSteps(numBlocks, 0.1, 0.8);
971
972
973
974
975

    // -------------------------------------------------------------------
    // LOAD THE DATA
    //--------------------------------------------------------------------
    PRAGMA_OMP( parallel for schedule(dynamic, 1) if (parallelProcessing) )
976
    for (int blockNum = 0; blockNum < int(numBlocks); blockNum++) {
977
978
979
980
981
      PARALLEL_START_INTERUPT_REGION

      // Find the workspace for this particular thread
      EventWorkspace_sptr ws;
      size_t threadNum = 0;
982
      if (parallelProcessing) {
983
984
        threadNum = PARALLEL_THREAD_NUMBER;
        ws = partWorkspaces[threadNum];
985
      } else
986
987
988
        ws = workspace;

      // Get the buffer (for this thread)
989
      DasEvent *event_buffer = buffers[threadNum];
990
991

      // Get the speeding-up array of vector<tofEvent> where index = detid.
992
      EventVector_pt *theseEventVectors = eventVectors[threadNum];
993
994
995
996
997

      // Where to start in the file?
      size_t fileOffset = first_event + (loadBlockSize * blockNum);
      // May need to reduce size of last (or only) block
      size_t current_event_buffer_size =
998
999
1000
          (blockNum == int(numBlocks - 1))
              ? (m_maxNumEvents - (numBlocks - 1) * loadBlockSize)
              : loadBlockSize;
For faster browsing, not all history is shown. View entire blame