LoadEventPreNexus2.cpp

  if (event_indices.size() < num_pulses) {
    g_log.warning()
        << "Event_indices vector is smaller than the pulsetimes array.\n";
    numPulses = static_cast<int64_t>(event_indices.size());
  }

  // Local stastic parameters
  size_t local_num_error_events = 0;
  size_t local_num_bad_events = 0;
  size_t local_num_wrongdetid_events = 0;
  size_t local_num_ignored_events = 0;
  size_t local_num_good_events = 0;
  double local_shortest_tof =
      static_cast<double>(MAX_TOF_UINT32) * TOF_CONVERSION;
  double local_longest_tof = 0.;

  // Storages
  std::map<PixelType, size_t> local_pidindexmap;
  std::vector<std::vector<Kernel::DateAndTime>> local_pulsetimes;
  std::vector<std::vector<double>> local_tofs;
  std::set<PixelType> local_wrongdetids;

  // process the individual events
  std::stringstream dbss;
  // size_t numwrongpid = 0;
  for (size_t i = 0; i < current_event_buffer_size; i++) {
    DasEvent &temp = *(event_buffer + i);
    PixelType pid = temp.pid;
    bool iswrongdetid = false;

    if (dbprint && i < m_dbOpNumEvents)
      dbss << i << " \t" << temp.tof << " \t" << temp.pid << "\n";

    // Filter out bad event
    if ((pid & ERROR_PID) == ERROR_PID) {
      local_num_error_events++;
      local_num_bad_events++;
      continue;
    }

    // Covert the pixel ID from DAS pixel to our pixel ID
    // downstream monitor pixel for SNAP
    if (pid == 1073741843)
      pid = 1179648;
    else if (this->using_mapping_file) {
      PixelType unmapped_pid = pid % this->numpixel;
      pid = this->pixelmap[unmapped_pid];
    }

    // Wrong pixel IDs
    if (pid > static_cast<PixelType>(detid_max)) {
      iswrongdetid = true;

      local_num_error_events++;
      local_num_wrongdetid_events++;
      local_wrongdetids.insert(pid);
    }

    // Now check if this pid we want to load.
    if (loadOnlySomeSpectra && !iswrongdetid) {
      std::map<int64_t, bool>::iterator it;
      it = spectraLoadMap.find(pid);
      if (it == spectraLoadMap.end()) {
        // Pixel ID was not found, so the event is being ignored.
        local_num_ignored_events++;
        continue;
      }
    }

    // Upon this point, only 'good' events are left to work on

    // Pulse: Find the pulse time for this event index
    if (pulse_i < numPulses - 1) {
      // This is the total offset into the file
      size_t total_i = i + fileOffset;
      // 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]))) {
        pulse_i++;
        if (pulse_i >= (numPulses - 1))
          break;
      }

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

    // TOF
    double tof = static_cast<double>(temp.tof) * TOF_CONVERSION;

    if (!iswrongdetid) {
      // Regular event that is belonged to a defined detector
      // Find the overall max/min tof
      if (tof < local_shortest_tof)
        local_shortest_tof = tof;
      if (tof > local_longest_tof)
        local_longest_tof = tof;

// This is equivalent to
// workspace->getEventList(this->pixel_to_wkspindex[pid]).addEventQuickly(event);
// But should be faster as a bunch of these calls were cached.
#if defined(__GNUC__) && !(defined(__INTEL_COMPILER)) && !(defined(__clang__))
      // This avoids a copy constructor call but is only available with GCC
      // (requires variadic templates)
      arrayOfVectors[pid]->emplace_back(tof, pulsetime);
#else
      arrayOfVectors[pid]->push_back(TofEvent(tof, pulsetime));
#endif

      ++local_num_good_events;
    } else {
      // Special events/Wrong detector id
      // i.  get/add index of the entry in map
      std::map<PixelType, size_t>::iterator it;
      it = local_pidindexmap.find(pid);
      size_t theindex = 0;
      if (it == local_pidindexmap.end()) {
        // Initialize it!
        size_t newindex = local_pulsetimes.size();
        local_pidindexmap[pid] = newindex;

        std::vector<Kernel::DateAndTime> tempvectime;
        std::vector<double> temptofs;
        local_pulsetimes.push_back(tempvectime);
        local_tofs.push_back(temptofs);

        theindex = newindex;

        // ++ numwrongpid;

        g_log.debug() << "Find New Wrong Pixel ID = " << pid << "\n";
      } else {
        // existing
        theindex = it->second;
      }

      // ii. calculate and add absolute time
      // int64_t abstime = (pulsetime.totalNanoseconds()+int64_t(tof*1000));
      local_pulsetimes[theindex].push_back(pulsetime);
      local_tofs[theindex].push_back(tof);

    } // END-IF-ELSE: On Event's Pixel's Nature

  } // ENDFOR each event

  if (dbprint)
    g_log.information(dbss.str());

  // Update local statistics to their global counterparts
  PARALLEL_CRITICAL(LoadEventPreNexus2_global_statistics) {
    this->num_good_events += local_num_good_events;
    this->num_ignored_events += local_num_ignored_events;
    this->num_error_events += local_num_error_events;

    this->num_bad_events += local_num_bad_events;
    this->num_wrongdetid_events += local_num_wrongdetid_events;

    std::set<PixelType>::iterator it;
    for (it = local_wrongdetids.begin(); it != local_wrongdetids.end(); ++it) {
      PixelType tmpid = *it;
      this->wrongdetids.insert(*it);

      // Create class map entry if not there
      size_t mindex = 0;
      auto git = this->wrongdetidmap.find(tmpid);
      if (git == this->wrongdetidmap.end()) {
        // create entry
        size_t newindex = this->wrongdetid_pulsetimes.size();
        this->wrongdetidmap[tmpid] = newindex;

        std::vector<Kernel::DateAndTime> temppulsetimes;
        std::vector<double> temptofs;
        this->wrongdetid_pulsetimes.push_back(temppulsetimes);
        this->wrongdetid_tofs.push_back(temptofs);

        mindex = newindex;
      } else {
        mindex = git->second;
      }

      // 2. Find local
      auto lit = local_pidindexmap.find(tmpid);
      size_t localindex = lit->second;

      for (size_t iv = 0; iv < local_pulsetimes[localindex].size(); iv++) {
        this->wrongdetid_pulsetimes[mindex].push_back(
            local_pulsetimes[localindex][iv]);
        this->wrongdetid_tofs[mindex].push_back(local_tofs[localindex][iv]);
      }
      // std::sort(this->wrongdetid_abstimes[mindex].begin(),
      // this->wrongdetid_abstimes[mindex].end());
    }

    if (local_shortest_tof < shortest_tof)
      shortest_tof = local_shortest_tof;
    if (local_longest_tof > longest_tof)
      longest_tof = local_longest_tof;
  } // END_CRITICAL

  return;
}

//----------------------------------------------------------------------------------------------
/** Comparator for sorting dasevent lists
  */
bool vzintermediatePixelIDComp(IntermediateEvent x, IntermediateEvent y) {
  return (x.pid < y.pid);
}

//-----------------------------------------------------------------------------
/**
  * 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)
  *
  * @param workspace :: Event workspace to set the proton charge on
  */
void LoadEventPreNexus2::setProtonCharge(
    DataObjects::EventWorkspace_sptr &workspace) {
  if (this->proton_charge.empty()) // nothing to do
    return;

  Run &run = workspace->mutableRun();

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

  // Add the time and associated charge to the log
  log->addValues(this->pulsetimes, this->proton_charge);

  /// TODO set the units for the log
  run.addLogData(log);
  // Force re-integration
  run.integrateProtonCharge();
  double integ = run.getProtonCharge();

  g_log.information() << "Total proton charge of " << integ
                      << " microAmp*hours found by integrating.\n";

  return;
}

//-----------------------------------------------------------------------------
/** Load a pixel mapping file
  * @param filename :: Path to file.
  */
void LoadEventPreNexus2::loadPixelMap(const std::string &filename) {
  this->using_mapping_file = false;
  this->pixelmap.clear();

  // check that there is a mapping file
  if (filename.empty()) {
    this->g_log.information("NOT using a mapping file");
    return;
  }

  // actually deal with the file
  this->g_log.debug("Using mapping file \"" + filename + "\"");

  // Open the file; will throw if there is any problem
  BinaryFile<PixelType> pixelmapFile(filename);
  PixelType max_pid = static_cast<PixelType>(pixelmapFile.getNumElements());
  // Load all the data
  pixelmapFile.loadAllInto(this->pixelmap);

  // Check for funky file
  if (std::find_if(pixelmap.begin(), pixelmap.end(),
                   std::bind2nd(std::greater<PixelType>(), max_pid)) !=
      pixelmap.end()) {
    this->g_log.warning("Pixel id in mapping file was out of bounds. Loading "
                        "without mapping file");
    this->numpixel = 0;
    this->pixelmap.clear();
    this->using_mapping_file = false;
    return;
  }

  // If we got here, the mapping file was loaded correctly and we'll use it
  this->using_mapping_file = true;
  // Let's assume that the # of pixels in the instrument matches the mapping
  // file length.
  this->numpixel = static_cast<uint32_t>(pixelmapFile.getNumElements());
}

//-----------------------------------------------------------------------------
/** Open an event file
  * @param filename :: file to open.
  */
void LoadEventPreNexus2::openEventFile(const std::string &filename) {
  // Open the file
  eventfile = new BinaryFile<DasEvent>(filename);
  num_events = eventfile->getNumElements();
  g_log.debug() << "File contains " << num_events << " event records.\n";

  // Check if we are only loading part of the event file
  const int chunk = getProperty("ChunkNumber");
  if (isEmpty(chunk)) // We are loading the whole file
  {
    first_event = 0;
    max_events = num_events;
  } else // We are loading part - work out the event number range
  {
    const int totalChunks = getProperty("TotalChunks");
    max_events = num_events / totalChunks;
    first_event = (chunk - 1) * max_events;
    // Need to add any remainder to the final chunk
    if (chunk == totalChunks)
      max_events += num_events % totalChunks;
  }

  g_log.information() << "Reading " << max_events << " event records\n";

  return;
}

//-----------------------------------------------------------------------------
/** Read a pulse ID file
 * @param filename :: file to load.
 * @param throwError :: Flag to trigger error throwing instead of just logging
 */
void LoadEventPreNexus2::readPulseidFile(const std::string &filename,
                                         const bool throwError) {
  this->proton_charge_tot = 0.;
  this->num_pulses = 0;
  this->pulsetimesincreasing = true;

  // jump out early if there isn't a filename
  if (filename.empty()) {
    this->g_log.information("NOT using a pulseid file");
    return;
  }

  std::vector<Pulse> *pulses;

  // set up for reading
  // Open the file; will throw if there is any problem
  try {
    BinaryFile<Pulse> pulseFile(filename);

    // Get the # of pulse
    this->num_pulses = pulseFile.getNumElements();
    this->g_log.information() << "Using pulseid file \"" << filename
                              << "\", with " << num_pulses << " pulses.\n";

    // Load all the data
    pulses = pulseFile.loadAll();
  } catch (runtime_error &e) {
    if (throwError) {
      throw;
    } else {
      this->g_log.information()
          << "Encountered error in pulseidfile (ignoring file): " << e.what()
          << "\n";
      return;
    }
  }

  if (num_pulses > 0) {
    double temp;
    DateAndTime lastPulseDateTime(0, 0);
    this->pulsetimes.reserve(num_pulses);
    for (size_t i = 0; i < num_pulses; i++) {
      Pulse &it = (*pulses)[i];
      DateAndTime pulseDateTime(static_cast<int64_t>(it.seconds), static_cast<int64_t>(it.nanoseconds));
      this->pulsetimes.push_back(pulseDateTime);
      this->event_indices.push_back(it.event_index);

      if (pulseDateTime < lastPulseDateTime)
        this->pulsetimesincreasing = false;
      else
        lastPulseDateTime = pulseDateTime;

      temp = it.pCurrent;
      this->proton_charge.push_back(temp);
      if (temp < 0.)
        this->g_log.warning("Individual proton charge < 0 being ignored");
      else
        this->proton_charge_tot += temp;
    }
  }

  this->proton_charge_tot = this->proton_charge_tot * CURRENT_CONVERSION;

  if (m_dbOpNumPulses > 0) {
    std::stringstream dbss;
    for (size_t i = 0; i < m_dbOpNumPulses; ++i)
      dbss << "[Pulse] " << i << "\t " << event_indices[i] << "\t "
           << pulsetimes[i].totalNanoseconds() << "\n";
    g_log.information(dbss.str());
  }

  // Clear the vector
  delete pulses;
}

//----------------------------------------------------------------------------------------------
/** Process input properties for purpose of investigation
  */
void LoadEventPreNexus2::processInvestigationInputs() {
  m_dbOpBlockNumber = getProperty("DBOutputBlockNumber");
  if (isEmpty(m_dbOpBlockNumber)) {
    m_dbOutput = false;
    m_dbOpBlockNumber = 0;
  } else {
    m_dbOutput = true;

    int numdbevents = getProperty("DBNumberOutputEvents");
    m_dbOpNumEvents = static_cast<size_t>(numdbevents);
  }

  int dbnumpulses = getProperty("DBNumberOutputPulses");
  if (!isEmpty(dbnumpulses))
    m_dbOpNumPulses = static_cast<size_t>(dbnumpulses);
  else
    m_dbOpNumPulses = 0;

  return;
}

} // namespace DataHandling
} // namespace Mantid