LoadEventPreNexus2.cpp

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

  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.;

  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
  for (size_t i = 0; i < current_event_buffer_size; i++)
  {
    DasEvent & temp = *(event_buffer + i);
    PixelType pid = temp.pid;
    bool iswrongdetid = false;

    if ((pid & ERROR_PID) == ERROR_PID) // marked as bad
    {
      local_num_error_events++;
      local_num_bad_events ++;
      continue;
    }

    //Covert the pixel ID from DAS pixel to our pixel ID
    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;
      }
    }

    // work with the good guys

    //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;
      }

      //if (pulsetimes[pulse_i] != pulsetime)    std::cout << pulse_i << " at " << pulsetimes[pulse_i] << "\n";

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

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

    if (!iswrongdetid){
      //Find the overall max/min tof
      if (tof < local_shortest_tof)
        local_shortest_tof = tof;
      if (tof > local_longest_tof)
        local_longest_tof = tof;

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

      // 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.
      arrayOfVectors[pid]->push_back(TofEvent(tof, pulsetime));

      // TODO work with period
      local_num_good_events++;

    } else {
      // b) 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;
      } 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);

    }

  } // ENDFOR each event

  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);

      // 1. Create class map entry if not there
      size_t mindex = 0;
      std::map<PixelType, size_t>::iterator 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
      std::map<PixelType, size_t>::iterator lit = local_pidindexmap.find(tmpid);
      size_t localindex= lit->second;

      // g_log.notice() << "Pixel " << tmpid << "  Global index = " << mindex << "  Local Index = " << localindex << std::endl;

      // 3. Sort and merge
      /* Redo
      std::sort(local_abstimes[localindex].begin(), local_abstimes[localindex].end());
      for (size_t iv = 0; iv < local_abstimes[localindex].size(); iv ++){
        this->wrongdetid_abstimes[mindex].push_back(local_abstimes[localindex][iv]);
      }
      */
      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;
  }
}

//-----------------------------------------------------------------------------
/// 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);
  double integ = run.integrateProtonCharge();
  //run.setProtonCharge(this->proton_charge_tot); //This is now redundant
  this->g_log.information() << "Total proton charge of " << integ << " microAmp*hours found by integrating.\n";

}

//-----------------------------------------------------------------------------
/** 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";
}

//-----------------------------------------------------------------------------
/** 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;
    }
  }

  double temp;

  if (num_pulses > 0)
  {
    DateAndTime lastPulseDateTime(0, 0);
    this->pulsetimes.reserve(num_pulses);
    for (size_t i=0; i < num_pulses; i++)
    {
      Pulse & it = (*pulses)[i];
      DateAndTime pulseDateTime( (int64_t) it.seconds, (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;

  //Clear the vector
  delete pulses;

}

} // namespace DataHandling
} // namespace Mantid