vDataFile.cc

Go to the documentation of this file.

/*
  Implementation of class vDataFile

  - separators: X
  - doxygen: X
  - default initialization: X
  - CASTOR_DEBUG: X
  - CASTOR_VERBOSE: X
*/

#include "vDataFile.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

vDataFile::vDataFile() 
{
  mp_ID = NULL;
  m_verbose = -1;
  
  // Variables related to the acquisition
  m2p_dataFile = NULL;
  m_headerFileName = "";
  m_dataFileName = "";
  m_scannerName = "";
  m_dataMode = MODE_UNKNOWN;
  m_dataType = TYPE_UNKNOWN;
  m_totalNbEvents = -1;
  m_dataType = -1;
  m_bedIndex = -1;
  m_startTimeInSec = 0.; 
  m_durationInSec = -1.;
  m_calibrationFactor = 1.;
  m_sizeEvent = -1;

  // Default POI (meaning we do not have any)
  mp_POIResolution[0] = -1.;
  mp_POIResolution[1] = -1.;
  mp_POIResolution[2] = -1.;
  mp_POIDirectionFlag[0] = false;
  mp_POIDirectionFlag[1] = false;
  mp_POIDirectionFlag[2] = false;
  m_POIInfoFlag = false;
  m_ignorePOIFlag = false;
  
  // Variable related to Buffer/Container arrays
  m2p_BufferEvent = NULL;
  mp_arrayEvents = NULL;
  m2p_bufferEventFromFile = NULL;
  m_mpi1stEvent = -1;
  m_mpiLastEvent = -1;
  m_mpiNbEvents = -1; 
  m_1stIdxArrayEvents = -1;
  m_lastIdxArrayEvents = -1;
  m_sizeArrayEvents = 0;
  m_percentageLoad = -1;
  m_requestBufferFilling = false;
  m_currentlyFillingBuffer = false;
  mp_currentlyReadingBuffer = NULL;
  mp_overBufferRange = NULL;
  mp_nbEventsReadFromFile = NULL;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

vDataFile::~vDataFile() 
{
  // Free array event buffer
  if (m_sizeArrayEvents > 0) if (mp_arrayEvents != NULL) delete mp_arrayEvents;
    
  // Free the (char**) m2p_bufferEventFromFile
  if (m2p_bufferEventFromFile)
  {
    for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
      if (m2p_bufferEventFromFile[th]) delete m2p_bufferEventFromFile[th];
    delete[] m2p_bufferEventFromFile;
  }
  // Free the (vEvent) buffer event
  if (m2p_BufferEvent)
  {
    for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
      if (m2p_BufferEvent[th]) delete m2p_BufferEvent[th];
    delete[] m2p_BufferEvent;
  }
  // Close datafiles and delete them
  if (m2p_dataFile)
  {
    for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
      if (m2p_dataFile[th]) m2p_dataFile[th]->close();
    delete m2p_dataFile;
  }
  // Destroy the boolean tabs that are used to organize and synchronize buffer reading
  if (mp_currentlyReadingBuffer) free(mp_currentlyReadingBuffer);
  if (mp_overBufferRange) free(mp_overBufferRange);
  // Destroy events reading counters, but print out before that
  if (mp_nbEventsReadFromFile)
  {
    // Verbose
    if (m_verbose>=4)
    {
      Cout("vDataFile::~vDataFile() -> Number of events directly read from file for each thread" << endl);
      int64_t total = 0;
      for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
      {
        Cout("  --> Thread " << th+1 << "  |  " << mp_nbEventsReadFromFile[th] << " events read from file" << endl);
        total += mp_nbEventsReadFromFile[th];
      }
      Cout("  --> Total number of events read from file: " << total << endl);
      Cout("  --> Number of events in the datafile: " << m_totalNbEvents << endl);
    }
    // Free memory
    free(mp_nbEventsReadFromFile);
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::ReadInfoInHeader()
{
  if (m_verbose>=3) Cout("vDataFile::ReadInfoInHeader() -> Read datafile header from '" << m_headerFileName << " ...'" << endl);
  
  // Read mandatory general fields in the header, check if errors (either mandatory tag not found or issue during data reading/conversion (==1) )
  if (ReadDataASCIIFile(m_headerFileName, "Number of events", &m_totalNbEvents, 1, KEYWORD_MANDATORY) )
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read number of events in the header data file !" << endl);
    return 1;
  }
  if (ReadDataASCIIFile(m_headerFileName, "Scanner name", &m_scannerName, 1, KEYWORD_MANDATORY) ) 
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read scanner name in the header data file !" << endl);
    return 1;
  }

  // Get data file name
  if (ReadDataASCIIFile(m_headerFileName, "Data filename", &m_dataFileName, 1, KEYWORD_MANDATORY) ) 
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data filename in the header data file !" << endl);
    return 1;
  }
  // If it is an absolute path (start with a OS_SEP) then let it as is, otherwise past the path of the header file name (because we suppose the two
  // files are side by side).
  // This is currently only unix compatible...
  if (m_dataFileName.substr(0,1)!=OS_SEP && m_headerFileName.find(OS_SEP)!=string::npos)
  {
    // Extract the path from the header file name
    size_t last_slash = m_headerFileName.find_last_of(OS_SEP);
    // Paste the path to the data file name
    m_dataFileName = m_headerFileName.substr(0,last_slash+1) + m_dataFileName;
  }

  // For data mode, multiple declarations are allowed, so we get the mode as a string and do some checks
  string data_mode = "";
  if (ReadDataASCIIFile(m_headerFileName, "Data mode", &data_mode, 1, KEYWORD_MANDATORY) ) 
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data mode in the header data file !" << endl);
    return 1;
  }
  if ( data_mode=="LIST" || data_mode=="LISTMODE" || data_mode=="LIST-MODE" ||
       data_mode=="list" || data_mode=="listmode" || data_mode=="list-mode" ||  data_mode=="0" ) m_dataMode = MODE_LIST;
  else if ( data_mode=="HISTOGRAM" || data_mode=="histogram" || data_mode=="Histogram" ||
            data_mode=="HISTO" || data_mode=="histo" || data_mode=="Histo" || data_mode=="1" ) m_dataMode = MODE_HISTOGRAM;
  else if ( data_mode=="NORMALIZATION" || data_mode=="normalization" || data_mode=="Normalization" ||
            data_mode=="NORM" || data_mode=="norm" || data_mode=="Norm" || data_mode=="2" ) m_dataMode = MODE_NORMALIZATION;
  else
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Unknown data mode '" << data_mode << "' found in header file !" << endl);
    return 1;
  }

  // For data type, multiple declarations are allowed, so we get the mode as a string and do some checks
  string data_type = "";
  if (ReadDataASCIIFile(m_headerFileName, "Data type", &data_type, 1, KEYWORD_MANDATORY) ) 
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data type in the header data file !" << endl);
    return 1;
  }
  if ( data_type=="PET" || data_type=="pet" || data_type=="0" ) m_dataType = TYPE_PET;
  else if ( data_type=="SPECT" || data_type=="spect" || data_type=="1" ) m_dataType = TYPE_SPECT;
  else if ( data_type=="TRANSMISSION" || data_type=="transmission" || data_type=="Transmission" ||
            data_type=="trans" || data_type=="TRANS" || data_type=="2" ) m_dataType = TYPE_TRANSMISSION;
  else
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Unknown data type '" << data_type << "' found in header file !" << endl);
    return 1;
  }

  // Get start time and duration, and then set the acquisition timing (do not do this for MODE_NORMALIZATION)
  if (m_dataMode!=MODE_NORMALIZATION)
  {
    if (ReadDataASCIIFile(m_headerFileName, "Start time (s)", &m_startTimeInSec, 1, KEYWORD_MANDATORY) ) 
    {
      Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read acquisition start time in the header data file !" << endl);
      return 1;
    }
    if (ReadDataASCIIFile(m_headerFileName, "Duration (s)", &m_durationInSec, 1, KEYWORD_MANDATORY) ) 
    {
      Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read acquisition stop time in the header data file !" << endl);
      return 1;
    }
    if (mp_ID->SetAcquisitionTime(m_bedIndex, m_startTimeInSec, m_durationInSec))
    {
      Cerr("***** vDataFile::ReadInfoInHeader() -> Error while passing acquisition start and stop time to the ImageDimensionsAndQuantification !" << endl);
      return 1;
    }
  }

  // Read optional fields in the header, check if errors (issue during data reading/conversion (==1) )
  if (ReadDataASCIIFile(m_headerFileName, "Calibration factor", &m_calibrationFactor, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "POI capability", mp_POIResolution, 3, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "POI correction flag", &m_POIInfoFlag, 3, KEYWORD_OPTIONAL)==1 )
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while reading optional field in the header data file !" << endl);
    return 1;
  }
  
  // Read fields specific to the modality (call to the related function implemented in child classes)
  if (ReadSpecificInfoInHeader() )
  {
    Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read modality-specific informations from the header data file !" << endl);
    return 1;
  }

  // Give the calibration factor to the oImageDimensionsAndQuantification that manages the quantification factors
  if (m_dataMode!=MODE_NORMALIZATION && mp_ID->SetCalibrationFactor(m_bedIndex, m_calibrationFactor))
  {
    Cerr("***** vDataFile::ReadSpecificInfoInHeader() -> A problem occured while setting the calibration factor to oImageDimensionsAndQuantification !" << endl);
    return 1;
  }

  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::CheckParameters()
{
  if(m_verbose >=3) Cout("vDataFile::CheckParameters() ..." << endl);

  if (mp_ID == NULL)
  {
    Cerr("***** vDataFile::CheckParameters() -> Error : ImageDimensionsAndQuantification object not initialized !" << endl);
    return 1;
  }
  if (m_headerFileName == "")
  {
    Cerr("***** vDataFile::CheckParameters() -> Error : string containing path to header file not initialized !" << endl);
    return 1;
  }  
  if (m_dataFileName == "")
  {
    Cerr("***** vDataFile::CheckParameters() -> Error : string containing path to raw data file not initialized !" << endl);
    return 1;
  }
  if (m_percentageLoad<0 || m_percentageLoad>100)
  {
    Cerr("***** vDataFile::CheckParameters() -> Percentage load of the data file incorrectly set !" << endl);
    return 1;
  }
  if (m_totalNbEvents<0)
  {
    Cerr("***** vDataFile::CheckParameters() -> Number of events incorrectly initialized !" << endl);
    return 1;
  }
  if (m_dataMode!=MODE_LIST && m_dataMode!=MODE_HISTOGRAM && m_dataMode!=MODE_NORMALIZATION)
  {
    Cerr("***** vDataFile::CheckParameters() -> Data mode incorrectly initialized !" << endl);
    return 1;
  }
  if (m_dataMode!= MODE_NORMALIZATION && m_durationInSec<0)
  {
    Cerr("***** vDataFile::CheckParameters() -> Acquisition duration (s) incorrectly initialized !" << endl);
    return 1;
  }
  if (m_dataType!=TYPE_PET && m_dataType!=TYPE_SPECT && m_dataType!=TYPE_TRANSMISSION)
  {
    Cerr("***** vDataFile::CheckParameters() -> Data type incorrectly initialized !" << endl);
    return 1;
  }
  if (m_bedIndex<0)
  {
    Cerr("***** vDataFile::CheckParameters() -> Bed position index incorrectly initialized !" << endl);
    return 1;
  }
  if (m_verbose<0)
  {
    Cerr("***** vDataFile::CheckParameters() -> Verbosity incorrectly initialized !" << endl);
    return 1;
  }
  
  // Call to the related function implemented in child classes
  if (CheckSpecificParameters())
  {
    Cerr("***** vDataFile::CheckParameters() -> Error while checking specific parameters !" << endl);
    return 1;
  }

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::InitializeFile()
{
  if (m_verbose>=3) Cout("vDataFile::InitializeFile() ..." << endl);
  
  // Instantiate as many file streams as threads
  m2p_dataFile = new fstream*[mp_ID->GetNbThreadsForProjection()];
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    m2p_dataFile[th] = new fstream( m_dataFileName.c_str(), ios::binary| ios::in );
    // If the datafile does not exist, notify to the user and exit program.
    if (!m2p_dataFile[th]->is_open()) 
    {
      Cerr("***** vDataFile::InitializeFile() -> Error reading the input data file at the path: '" << m_dataFileName.c_str() << "'" << endl);
      Cerr("                                     (provided in the data file header: " << m_headerFileName << ")" << endl);
      return 1;
    }
  }
  // Check file size consistency here (this is a pure virtual function implemented by children)
  if (CheckFileSizeConsistency())
  {
    Cerr("***** vDataFile::InitializeFile() -> A problem occured while checking file size consistency !" << endl);
    return 1;
  }
  // Initialization of the buffer related to direct data reading from file (during reconstruction)
  m2p_bufferEventFromFile = new char*[mp_ID->GetNbThreadsForProjection()];
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    m2p_bufferEventFromFile[th] = new char[m_sizeEvent];
  }
  // Initialization of the boolean tab used to know if any thread is currently reading into the file buffer
  mp_currentlyReadingBuffer = (bool*)malloc(mp_ID->GetNbThreadsForProjection()*sizeof(bool));
  for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) mp_currentlyReadingBuffer[th] = false;
  // Initialization of the boolean tab used to know if some threads are under the file buffer range
  mp_overBufferRange = (bool*)malloc(mp_ID->GetNbThreadsForProjection()*sizeof(bool));
  for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) mp_overBufferRange[th] = false;
  // Initialization of the counters for how many events are read directly from file
  mp_nbEventsReadFromFile = (int64_t*)malloc(mp_ID->GetNbThreadsForProjection()*sizeof(int64_t));
  for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) mp_nbEventsReadFromFile[th] = 0;

  // ------------------ Compute here the piece of file that each MPI instance manages --------------------

  // 1. Compute the number of events that each instance has to manage
  int64_t instance_size = m_totalNbEvents / mp_ID->GetMPISize();
  // All instances manage an equal part of the file but the last one which also manages the few last events
  if (mp_ID->GetMPIRank()!=mp_ID->GetMPISize()-1) m_mpiNbEvents = instance_size;
  else m_mpiNbEvents = instance_size + (m_totalNbEvents - instance_size*mp_ID->GetMPISize());
  // 2. Compute the first event managed by the instance
  m_mpi1stEvent = mp_ID->GetMPIRank() * instance_size;
  // 3. Compute the last event managed by the instance
  m_mpiLastEvent = m_mpi1stEvent + m_mpiNbEvents - 1;
    
  // Instanciation of the array of preloaded events (if enabled)
  if (m_percentageLoad>0)
  {
    if (m_verbose>=3)
    {
      Cout("  --> Allocating memory for the datafile" << endl);
      Cout("  --> Number of events in memory: " << m_mpiNbEvents << endl);
    }
    // Computation of the exact size of the buffer of preloaded events from the provided percentage load
    m_sizeArrayEvents = (((int64_t)m_percentageLoad) * m_mpiNbEvents) / 100;
    // Compute the number of bytes
    int64_t nb_bytes_to_load = m_sizeArrayEvents*m_sizeEvent;
    // Verbose to say that we proceed to the first reading of the datafile
    if (m_verbose>=3) Cout("  --> Load of " << m_percentageLoad << "% of the datafile (" << nb_bytes_to_load << " bytes)" << endl);
    // Instantiate the array of events in which the file content is directly copied
    mp_arrayEvents = new char [nb_bytes_to_load];
  }
  else m_sizeArrayEvents = 0;

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

void vDataFile::ResetBufferRange()
{
  // Reset the range indices only if the percentage is strictly between 0 and 100.
  // For the case 0, it is simply a security.
  // For the case 100, it would force the buffer to be re-read again...
  // The purpose of this function is to be called before a loop over the datafile, so that all
  // threads will have to deal with event indices always inside or above the range (never under).
  // This is important to avoid that the thread 0 filling the buffer is ahead of other threads, leading to the other threads
  // not benefiting from the file buffer. There is a mechanism in the GetEvent() function so that the thread 0 waits for the
  // other threads to be over the buffer range before loading a new portion of the data file into the buffer.
  if (m_percentageLoad>0 && m_percentageLoad<100)
  {
    m_1stIdxArrayEvents = -1;
    m_lastIdxArrayEvents = -1;
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::FillBuffer(int64_t a_eventIndex)
{
  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::FillBuffer() ..." << endl);
  #endif

  // Important note: m_mpiLastEvent/m_lastIdxArrayEvent are the last event indices INCLUDED in the instance/buffer

  // Note: This function is called from the GetEvent() function inside the iterations, but only by thread #0
  int thread_0 = 0;

  // Seek to the appropriate position in the input data file
  m2p_dataFile[thread_0]->seekg(m_sizeEvent*a_eventIndex);

  // Compute the last event index that can be loaded into the buffer (included; not excluded)
  int64_t last_included_event_that_can_be_loaded = a_eventIndex + m_sizeArrayEvents - 1;

  // This is the number of events that will be loaded
  int64_t nb_events_to_load = 0;
  // If the number of remaining events to be loaded is less than the buffer capacity, we only load until m_mpiLastEvent is reached
  if (last_included_event_that_can_be_loaded > m_mpiLastEvent) nb_events_to_load = m_mpiLastEvent + 1 - a_eventIndex;
  // Otherwise we fully load the buffer
  else nb_events_to_load = m_sizeArrayEvents;

  // Read the data
  m2p_dataFile[thread_0]->read(mp_arrayEvents, nb_events_to_load*m_sizeEvent);

  // Todo: histogram/SPECT only!!!!!!!!!!!!!!
  //Shuffle( nb_events_to_load );

  // Check if the reading was successfull
  if (!(*m2p_dataFile[thread_0]))
  {
    Cerr("***** vDataFile::FillBuffer() -> Failed to read " << nb_events_to_load << " events in datafile (only read " << m2p_dataFile[thread_0]->gcount() << " events) !" << endl);
    return 1;
  }

  // Update the buffer first and last indices
  m_1stIdxArrayEvents  = a_eventIndex;
  m_lastIdxArrayEvents = m_1stIdxArrayEvents + nb_events_to_load - 1;

  // Verbose
  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::FillBuffer() completed" << endl);
  #endif

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::Shuffle( int64_t nb_events_to_load )
{
  // Buffer storing ordered indices from ( 0 ) to ( nb_events_to_load - 1 )
  int64_t *rndmIdx = new int64_t[ nb_events_to_load ];
  std::iota( rndmIdx, rndmIdx + nb_events_to_load, 0 );

  // Initializing random
  std::random_device rd;
  std::mt19937 rndm( rd() );
  rndm.seed( 1100001001 );

  // Shuffling the buffer of indices
  std::shuffle( rndmIdx, rndmIdx + nb_events_to_load, rndm );

  // Creating a tmp buffer
  char *mp_arrayEvents_tmp = new char[ nb_events_to_load*m_sizeEvent ];

  // Copy sorted buffer to shuffled buffer
  for( int64_t i = 0; i < nb_events_to_load; ++i )
  {
    for( int64_t j = 0; j < m_sizeEvent; ++j )
    {
      mp_arrayEvents_tmp[ i * m_sizeEvent + j ] = mp_arrayEvents[ rndmIdx[ i ] * m_sizeEvent + j ];
    }
  }

  // Freeing memory
  delete[] rndmIdx;
  delete[] mp_arrayEvents;

  mp_arrayEvents = mp_arrayEvents_tmp;

  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

vEvent* vDataFile::GetEventWithAscendingOrderAssumption(int64_t a_eventIndex, int64_t a_eventIndexLimit, int a_th)
{
  // With the assumption that this function is called from a loop over event indices in an ascending order, the implementation
  // is as follows:
  //  - The ResetBufferRange() function MUST be called before launching the loop over event indices
  //  - Thread 0 is a particular case: it is the only one able to fill the buffer, and to do that it waits for the other threads to be above the current range.
  //  - The other threads read directly from the file if the first thread is currently filling the buffer, or if below or above the current range. In the latter
  //    case, they set a boolean flag saying they are above so that the thread 0 can know it. Otherwise they read from the buffer.
  // See the comments at the beginning of the GetEventWithoutOrderAssumption() function for details about motivations, etc.

  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::GetEventWithAscendingOrderAssumption() ..." << endl);
  #endif

  // Particular case for an empty buffer (-load 0): read the event directly from the file, whatever the thread number
  if (m_sizeArrayEvents==0) return GetEventFromFile(a_eventIndex, a_th);

  // Special treatment for the first thread (it is the only one allowed to refill the buffer)
  if (a_th==0)
  {
    // If the event index is below the current indices range loaded in memory, it is not normal as only this thread refills the buffer
    if (a_eventIndex<m_1stIdxArrayEvents)
    {
      Cerr("***** vDataFile::GetEventWithAscendingOrderAssumption() -> The requested event index is below the current range of indices loaded in memory !" << endl);
      Cerr("  --> Note that this situation should never happen if this function is being called within an ascending ordered loop on event indices, AND if" << endl);
      Cerr("      the vDataFile::ResetBufferRange() function has been called before performing the loop. So if this is really the case, then please report" << endl);
      Cerr("      the bug." << endl);
      return NULL;
    }
    // Else if the event is above the current indices range, then this first thread has to refill the buffer
    else if (a_eventIndex>m_lastIdxArrayEvents)
    {
      // Here, the thread has to wait for all others to be over the current buffer range
      bool all_threads_over_range = true;
      for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++)
        all_threads_over_range = all_threads_over_range && mp_overBufferRange[th];
      while (!all_threads_over_range)
      {
        all_threads_over_range = true;
        for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++)
          all_threads_over_range = all_threads_over_range && mp_overBufferRange[th];
      }
      // Now we are good, so we notify that the first thread is currently filling the buffer
      m_currentlyFillingBuffer = true;
      // Now the first thread can safely fill the buffer from the current event index, because all threads are above the range and they know the first
      // thread will refill the buffer
      if (FillBuffer(a_eventIndex))
      {
        Cerr("***** vDataFile::GetEventWithAscendingOrderAssumption() -> An error occured while filling the buffer from datafile !" << endl);
        return NULL;
      }
      // Say that the first thread is done filling the buffer
      m_currentlyFillingBuffer = false;
    }
    // Get address of the event to recover from the array buffer
    char* event_address_in_array = mp_arrayEvents + (a_eventIndex-m_1stIdxArrayEvents)*m_sizeEvent;
    // And read from the buffer
    return GetEventFromBuffer(event_address_in_array, a_th);
  }
  // For the other threads
  else
  {
    // The event pointer that will be returned
    vEvent* p_event;
    // We first check if the thread 0 is currently filling the buffer; in that case then read from file
    if (m_currentlyFillingBuffer)
    {
      // In this case, we still say that the thread is over the buffer range. In case this thread finishes its part of the loop
      // while the first thread is filling the buffer, the first thread won't be blocked into waiting for all threads to be over
      // the range.
      mp_nbEventsReadFromFile[a_th]++;
      p_event = GetEventFromFile(a_eventIndex, a_th);
    }
    // Else if the event index is below the current indices range
    else if (a_eventIndex<m_1stIdxArrayEvents)
    {
      // This situation may still occur in some cases: right after a buffer filling and if the number of subsets is high. This
      // is because the new portion of the data file loaded by the thread 0 is not strictly contiguous to the previous portion
      // due to subsets and multi-threading (in case of a block thread sequence higher than 1 or a high number of threads). In
      // this case we read directly from the file.
      mp_overBufferRange[a_th] = false;
      mp_nbEventsReadFromFile[a_th]++;
      p_event = GetEventFromFile(a_eventIndex, a_th);
    }
    // Else if the event index is above the current indices range
    else if (a_eventIndex>m_lastIdxArrayEvents)
    {
      // We say that this thread is above the current range, so that the thread 0 knows that
      mp_overBufferRange[a_th] = true;
      // And we read directly from the file
      mp_nbEventsReadFromFile[a_th]++;
      p_event = GetEventFromFile(a_eventIndex, a_th);
    }
    // Else, we can safely read from the buffer
    else
    {
      // We check again here that the first thread is not reading into the buffer in case it started between now and the first check.
      // This allows to catch some potential bad behaviors in extrem conditions (a lot of threads with a very small load percentage).
      if (m_currentlyFillingBuffer)
      {
        mp_nbEventsReadFromFile[a_th]++;
        p_event = GetEventFromFile(a_eventIndex, a_th);
      }
      else
      {
        // Say that we are not above range
        mp_overBufferRange[a_th] = false;
        // Compute the adress into the buffer
        char* event_address_in_array = mp_arrayEvents + (a_eventIndex-m_1stIdxArrayEvents)*m_sizeEvent;
        // Get the event
        p_event = GetEventFromBuffer(event_address_in_array, a_th);
      }
    }
    // If the current index is over the provided limit index, we say that this thread is over the buffer range to avoid thread 0 to be stuck
    if (a_eventIndex >= a_eventIndexLimit) mp_overBufferRange[a_th] = true;
    // Return the event
    return p_event;
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

vEvent* vDataFile::GetEventWithoutOrderAssumption(int64_t a_eventIndex, int a_th)
{
  // Without any assumption about the order of event indices, the implementation of this function is quite
  // tricky when we want genericity and efficiency! We absolutely want that only one thread deals with the
  // buffer filling (the first thread) while others are not blocked during this process and can still read
  // events directly from the file. Under these restrictions, multiple strategies were tested, but the current
  // one is the most efficient and stable one. The most stressful situation for this function is when using
  // many threads with a small load percentage. The use of OpenMP locks, in many different ways has always
  // lead to rare occurences (1-2%) of segmentation fault (load percentage typically between 1 to 10 and 16
  // threads), which could not be explained but also not tolerated... The current implementation uses home-made
  // locks at multiple levels but which does not hamper the efficiency. It has been validated with intense testing
  // (1500 reconstructions of the PET benchmark with a number of threads randomly chosen between 50 to 120 and a load
  // percentage of 2%; note that if somebody has such a powerfull computer with such a small amount of RAM, it would
  // be quite a paradoxal situation, but nonetheless, we wanted to stress this function at most). The results are:
  // 1 segmentation fault and 2 runs with a few pixels variating up to 0.2% (basically, it is one event at one moment
  // that succeeds in passing between the test of m_requestBufferFilling and the affectation of mp_currentlyReadingBuffer,
  // so that it runs through the GetEventFromBuffer function while the first thread also succeed to pass between the
  // checks and is filling the buffer at the same time; plus in order to give a corrupted event, the buffer has to be
  // changed right before the thread read at this memory location; with such rare stuff in such rare conditions, we feel
  // this function is sufficiently robust with respect to the restrictions we assume at the beginning; still, if
  // someone has something better, we take it!).
  // Anyway, if you want to get events inside an ascending ordered loop over events, then use the GetEventWithAscendingOrderAssumption()
  // function instead. It is more efficient and has the same failure rate in the same chaotic conditions (a lot of threads with a very
  // small load percentage. But do not forget to call the ResetBufferRange() function before the loop. However, note that this function
  // is not yet compatible with the use of MPI.

  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::GetEventWithoutOrderAssumption() ..." << endl);
  #endif

  // Particular case for an empty buffer (-load 0): read the event directly from the file, whatever the thread number
  if (m_sizeArrayEvents==0) return GetEventFromFile(a_eventIndex, a_th);

  // Special treatment for the first thread (it is the only one allowed to refill the buffer)
  if (a_th==0)
  {
    // If the event is outside the current buffer, then the thread refills the buffer
    if (a_eventIndex<m_1stIdxArrayEvents || a_eventIndex>m_lastIdxArrayEvents)
    {
      // Say that the first thread is requesting to fill the buffer
      m_requestBufferFilling = true;
      // Now this thread will wait for all other threads to declare that they are not currently reading into the buffer (otherwise
      // the reading of the event by the other threads may be corrupted). The other thread are now aware that the first thread is
      // waiting so they will not read from the buffer once finished their current reading.
      bool all_threads_not_reading_buffer = true;
      for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++)
        all_threads_not_reading_buffer = all_threads_not_reading_buffer && !mp_currentlyReadingBuffer[th];
      while (!all_threads_not_reading_buffer)
      {
        all_threads_not_reading_buffer = true;
        for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++)
          all_threads_not_reading_buffer = all_threads_not_reading_buffer && !mp_currentlyReadingBuffer[th];
      }
      // We are good, so say that the first thread is currently filling the buffer (this is needed below, see explanations)
      m_currentlyFillingBuffer = true;
      // Now the first thread can safely fill the buffer from the current event index, because as all threads know it was waiting for
      // filling the buffer, they know read directly from the file
      if (FillBuffer(a_eventIndex))
      {
        Cerr("***** vDataFile::GetEventWithoutOrderAssumption() -> An error occured while filling the buffer from datafile !" << endl);
        return NULL;
      }
      // Say that the first thread is done filling the buffer
      m_currentlyFillingBuffer = false;
      // No more request to fill the buffer
      m_requestBufferFilling = false;
    }
    // Get address of the event to recover from the array buffer
    char* event_address_in_array = mp_arrayEvents + (a_eventIndex-m_1stIdxArrayEvents)*m_sizeEvent;
    // And read from the buffer
    return GetEventFromBuffer(event_address_in_array, a_th);
  }
  // For the other threads
  else
  {
    // If out of buffer range, then read from file (we must check at first if the first thread is currently filling the buffer, otherwise
    // the range 1st and last indices may not reflect the actual range; then we cannot use here the m_requestBufferFilling flag instead,
    // because some compiler optimization might merge this if section with the next one, and we absolutely need the next one to be alone,
    // and right befre the first instruction of the else section setting the mp_currentlyReadingBuffer[a_th] flag to true.
    if (m_currentlyFillingBuffer || a_eventIndex<m_1stIdxArrayEvents || a_eventIndex>m_lastIdxArrayEvents)
    {
      mp_nbEventsReadFromFile[a_th]++;
      return GetEventFromFile(a_eventIndex, a_th);
    }
    // Else if request filling, then read from file (this check has to be separated from the one above, otherwise if there is no buffer
    // filling request, between the time this thread do the checks and set the currentlyReadingBuffer to true, the thread0 may have
    // enaugh time to start the filling of the buffer, creating erroneous lines, also may be due to some compiler optimizations reversing
    // some checks; at least, this version is stable)
    else if (m_requestBufferFilling)
    {
      mp_nbEventsReadFromFile[a_th]++;
      return GetEventFromFile(a_eventIndex, a_th);
    }
    // Else, we can safely read from the buffer
    else
    {
      // Say that this thread is reading into the buffer
      mp_currentlyReadingBuffer[a_th] = true;
      // Compute the adress into the buffer
      char* event_address_in_array = mp_arrayEvents + (a_eventIndex-m_1stIdxArrayEvents)*m_sizeEvent;
      // Get the event
      vEvent* p_event = GetEventFromBuffer(event_address_in_array, a_th);
      // Say that the thread has done reading into the buffer
      mp_currentlyReadingBuffer[a_th] = false;
      // Return the event
      return p_event;
    }
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

vEvent* vDataFile::GetEventFromFile(int64_t a_eventIndex, int a_th)
{
  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::GetEventFromFile() ..." << endl);
  #endif
  
  // Seek to the position and read the event
  m2p_dataFile[a_th]->seekg((int64_t)m_sizeEvent*a_eventIndex);
  m2p_dataFile[a_th]->read(m2p_bufferEventFromFile[a_th], m_sizeEvent);
  // Check the reading
  if (!(*m2p_dataFile[a_th]))
  {
    Cerr("***** vDataFile::GetEventFromFile() -> Failed to read one event in datafile !" << endl);
    return NULL;
  }
  // Pass the pointer to the event to the interpretor function (implemented by children)
  return GetEventFromBuffer(m2p_bufferEventFromFile[a_th], a_th);
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

void vDataFile::GetEventIndexStartAndStop( int64_t* ap_indexStart, int64_t* ap_indexStop, int a_subsetNum, int a_nbSubsets )
{
  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::GetEventIndexStartAndStop() ..." << endl);
  #endif
  // Basically here, the index start for a single MPI instance will be the current subset number.
  // If multiple instances are used, the whole datafile is split into equal-size concatenated pieces.
  // So for each instance, we just have to found the first index falling in its range (assuming that
  // the event index step is always equal to the number of subsets).

  // For the first machine, index start is the current subset number
  if (mp_ID->GetMPIRank()==0) *ap_indexStart = a_subsetNum;
  // For the other machines, we must search for the first index falling in their range belonging to this
  // subset (a subset being starting at a_subsetNum with a step equal to the number of subsets, a_nbSubsets)
  else
  {
    // Compute the modulo of the first index of this machine minus the subset number with respect to the number of subsets
    int64_t modulo = (m_mpi1stEvent-a_subsetNum) % a_nbSubsets;
    // If this modulo is null, then the index start is the first index
    if (modulo==0) *ap_indexStart = m_mpi1stEvent;
    // Otherwise, the index start is equal to the first index plus the number of subsets minus the modulo
    else *ap_indexStart = m_mpi1stEvent + (a_nbSubsets - modulo);
  }

  // For index stop, we simply get the last event of the MPI instance (+1 because the for loop is exclusive)
  *ap_indexStop = m_mpiLastEvent + 1;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int vDataFile::GetMaxRingDiff()
{
  Cerr("*****vDataFile::GetMaxRingDiff() ->This function is not implemented for the used system" << endl);
  Cerr("                                  (this error may be prompted if the present function is erroneously called for a SPECT system)" << endl);
  return -1;
}

/************************************************
* PROJECTION SCRIPT FUNCTIONS
************************************************/

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteData()
  \brief Write/Merge chunk of data in a general data file. 
  \todo check if some file manipulations are correctly done
  \todo adapt to the data loading/writing in RAM
  \return 0 if success, and positive value otherwise.
*/
int vDataFile::PROJ_WriteData()
{
  if (m_verbose >=2) Cout("vDataFile::PROJ_WriteData() ..." << endl);
  
  // Close all multithreaded datafiles before merging
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
    if (m2p_dataFile[th]) m2p_dataFile[th]->close();


  // If the output file doesn't exist yet, simply rename the first temporary file name using the ouput file name
  if (!ifstream(m_dataFileName))
  {
    // If only one projection is required (i.e one threads && no projection of frame or gates)...
    if(mp_ID->GetNbThreadsForProjection()==1  // No multithreading so no multiple tmp datafiles 
       && mp_ID->GetNbTimeFrames()*
       mp_ID->GetSensNbRespGates()*
       mp_ID->GetSensNbCardGates()  == 1)
       {
         // rename the first temporary file name to the global name
         string tmp_file_name = m_dataFileName + "_0";
    
         // ... then just rename the first temporary file name using the ouput file name
         rename(tmp_file_name.c_str(),m_dataFileName.c_str());
         // no need to concatenate, so we leave here.
         return 0 ;
       }
  }

  // Create the final output file which will concatenate the events inside the thread-specific data files
  ofstream merged_file(m_dataFileName.c_str(), ios::out | ios::binary | ios::app);

  // Concatenation : generate input file ("data_file") to read the buffer of the thread-specific data files and store the information in the final output file
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    // Build thread file name
    stringstream ss; ss << th;
    string file_name = m_dataFileName;
    file_name.append("_").append(ss.str());
    // Open it
    // Note SS: Maybe we can use the m2p_dataFile[th] here by just rewarding to the beginning of the file ?
    // Note TM: There were some issues involving the use of rdbuf and concatenation of the file in this case (ifstream were needed instead of fstream for some reasons)
    //          But we should have another look on how the projection data writing works with the implementation of the new analytical simulator.
    ifstream data_file(file_name.c_str(), ios::binary | ios::in);

    if (!data_file)
    {
      Cerr(endl);
      Cerr("***** vDataFile::PROJ_ConcatenateMthDatafile() -> Input temporary thread file '" << file_name << "' is missing or corrupted !" << endl);
      return 1;
    }

    // Concatenate it to the merged file
    merged_file << data_file.rdbuf();
    // Close file
    data_file.close();
    
    // Re-open datafiles (needed if projecting frame/gates, as the contents of the temporary datafile are copied to the main datafile after each frame/gate)
    m2p_dataFile[th]->open( file_name.c_str(), ios::binary | ios::out | ios::trunc);
  }

  // Close merged file
  merged_file.close();

  return 0;
}






// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_DeleteTmpDatafile()
  \brief  Delete temporary datafile used for multithreaded output writing if needed
  \return 0 if success, and positive value otherwise.
*/
int vDataFile::PROJ_DeleteTmpDatafile()
{
  if (m_verbose >=3) Cout("vDataFile::PROJ_DeleteTmpDatafile() ..." << endl);
  
  // Generate input file ("data_file") to read the buffer of the thread-specific data files and store the information in the final output file
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    // Build thread file name
    stringstream ss; ss << th;

    if (m2p_dataFile[th]) m2p_dataFile[th]->close();

    string file_name = m_dataFileName;
    file_name.append("_").append(ss.str());

    // Remove temporary file for data output writing (Projection script only)
    ifstream fcheck(file_name.c_str());
    if(fcheck.good())
    {
      fcheck.close();
      #ifdef _WIN32
      string dos_instruction = "del " + file_name;
      system(dos_instruction.c_str());
      #else
      remove(file_name.c_str());
      #endif
    }
  }
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_GenerateEvent(int idx_elt1, int idx_elt2, int a_th)
  \param idx_elt1 : first ID of the event
  \param idx_elt2 : second ID of the event
  \param a_th : index of the thread from which the function was called
  \brief  Generate a standard event and set up its ID
          Used by the projection, list-mode sensitivity generation, and datafile converter scripts
  \return the thread specific m2p_BufferEvent array containing the event
*/
vEvent* vDataFile::PROJ_GenerateEvent(int a_idxElt1, int a_idxElt2, int a_th)
{
  #ifdef CASTOR_VERBOSE
  if (m_verbose >=4) Cout("vDataFile::PROJ_GenerateEvent() ..." << endl);
  #endif
  
  // Only 1 line required for projection/sensitivity generation
  m2p_BufferEvent[a_th]->SetNbLines(1);  
  m2p_BufferEvent[a_th]->SetID1(0, a_idxElt1);
  m2p_BufferEvent[a_th]->SetID2(0, a_idxElt2);
  
  return m2p_BufferEvent[a_th];
}