iDataFilePET.cc

Go to the documentation of this file.

/*
  Implementation of class iDataFilePET

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

#include "iDataFilePET.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iDataFilePET constructor. 
         Initialize the member variables to their default values.
*/
iDataFilePET::iDataFilePET() : vDataFile() 
{
  // Read optional fields in the header
  m_dataType = TYPE_PET;
  m_maxNumberOfLinesPerEvent = 1;
  m_maxRingDiff = -1;
  m_isotope = "unknown";
  m_resolutionTOF = -1.;
  m_nbBinsTOF = -1;
  m_binSizeTOF = -1.;
  m_eventKindFlag = false;
  m_atnCorrectionFlag = false;
  m_ignoreAttnCorrectionFlag = false;
  m_normCorrectionFlag = false;
  m_ignoreNormCorrectionFlag = false;
  m_scatCorrectionFlag = false;
  m_ignoreScatCorrectionFlag = false;
  m_randCorrectionFlag = false;
  m_ignoreRandCorrectionFlag = false;
  m_TOFInfoFlag = false;
  m_ignoreTOFFlag = false;
}


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iDataFilePET destructor. 
*/
iDataFilePET::~iDataFilePET() {}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadSpecificInfoInHeader()
  \brief Read through the header file and gather specific PET information. 
  \return 0 is success, positive value otherwise
*/
int iDataFilePET::ReadSpecificInfoInHeader()
{ 
  if(m_verbose >=3) Cout("iDataFilePET::ReadSpecificInfoInHeader()... " << endl);
  
  // Read optional fields in the header, check if errors (issue during data reading/conversion (==1) )
  if (ReadDataASCIIFile(m_headerFileName, "Maximum number of lines per event", &m_maxNumberOfLinesPerEvent, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Maximum ring difference", &m_maxRingDiff, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "TOF resolution", &m_resolutionTOF, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "TOF correction flag", &m_TOFInfoFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Histo TOF bin", &m_nbBinsTOF, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Histo TOF size", &m_binSizeTOF, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Coincidence kind flag", &m_eventKindFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Attenuation correction flag", &m_atnCorrectionFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Normalization correction flag", &m_normCorrectionFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Scatter correction flag", &m_scatCorrectionFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Random correction flag", &m_randCorrectionFlag, 1, KEYWORD_OPTIONAL)==1 ||
      ReadDataASCIIFile(m_headerFileName, "Isotope", &m_isotope, 1, KEYWORD_OPTIONAL)==1 )
  {
    Cerr("***** iDataFilePET::ReadSpecificInfoInHeader() -> Error while reading optional field in the header data file '" << endl);
    return 1;
  }

  // Give the PET isotope to the oImageDimensionsAndQuantification that manages the quantification factors
  if (m_dataMode!=MODE_NORMALIZATION && mp_ID->SetPETIsotope(m_bedIndex, m_isotope))
  {
    Cerr("***** iDataFilePET::ReadSpecificInfoInHeader() -> A problem occured while setting the isotope to oImageDimensionsAndQuantification !" << endl);
    return 1;
  }

  return 0;
}


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ComputeSizeEvent()
  \brief Computation of the size of each event according to the mandatory/optional correction fields
  \todo Additionnal optionnal fields (physiological info (Phy field) & energy info)
  \return 0 is success, positive value otherwise
*/
int iDataFilePET::ComputeSizeEvent()
{
  // Verbose
  if (m_verbose>=3) Cout("iDataFilePET::ComputeSizeEvent() ..." << endl);

  // For MODE_LIST events
  if (m_dataMode == MODE_LIST) 
  {
    // Size of the mandatory element in a list-mode event:
    
    // Time + max_nb_lines*2*crystalID
    m_sizeEvent = sizeof(uint32_t) + m_maxNumberOfLinesPerEvent*2*sizeof(uint32_t);
    
    // Number of lines
    if (m_maxNumberOfLinesPerEvent>1) m_sizeEvent += sizeof(uint16_t);

    // Optional flags
    if (m_eventKindFlag) m_sizeEvent += sizeof(uint8_t);
    if (m_atnCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_scatCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_randCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_normCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_TOFInfoFlag) m_sizeEvent += sizeof(uint32_t);
    
    // POI availability in each direction
    for (int i=0 ; i<3; i++) if (mp_POIDirectionFlag[i]) m_sizeEvent += 2*sizeof(FLTNBDATA);
  }
  // For MODE_HISTOGRAM events
  else if (m_dataMode == MODE_HISTOGRAM)
  {
    // Size of the mandatory element in a histo event:
    
    // Time + nbBinsTOF*event_value + max_nb_line*2*crystalID
    m_sizeEvent = sizeof(uint32_t) + m_nbBinsTOF*sizeof(FLTNBDATA) + m_maxNumberOfLinesPerEvent*2*sizeof(uint32_t);
    
    // Number of lines
    if (m_maxNumberOfLinesPerEvent>1) m_sizeEvent += sizeof(uint16_t);
    
    // Optional flags
    if (m_atnCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_randCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_normCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_scatCorrectionFlag) m_sizeEvent += m_nbBinsTOF*sizeof(FLTNBDATA);
  }
  // For MODE_NORMALIZATION events
  else if (m_dataMode == MODE_NORMALIZATION)
  {
    // max_nb_lines*2*crystalID
    m_sizeEvent = m_maxNumberOfLinesPerEvent*2*sizeof(uint32_t);
    // Number of lines
    if (m_maxNumberOfLinesPerEvent>1) m_sizeEvent += sizeof(uint16_t);
    // Optional flag for normalization
    if (m_normCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    // Optional flag for attenuation
    if (m_atnCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
  }
  // Unknown event type
  else
  {
    Cerr("***** iDataFilePET::ComputeSizeEvent() -> Unknown event mode !" << endl);
    return 1;
  }

  if(m_verbose >=3) Cout("iDataFilePET::ComputeSizeEvent() -> Event size = " << m_sizeEvent << " bytes" << endl);

  if (m_sizeEvent <= 0) 
  {
    Cerr("***** iDataFilePET::ComputeSizeEvent() -> Error, the Event size in bytes should be >= 0 !" << endl;);
    return 1;
  }

  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PrepareDataFile()
  \brief Store different kind of information inside arrays (data relative to specific correction as well as basic raw data for the case data is loaded in RAM)
         Use the flag provided by the user to determine how the data has to be sorted (preloaded or read on the fly)
  \return 0 is success, positive value otherwise
*/
int iDataFilePET::PrepareDataFile()
{
  // Verbose
  if (m_verbose>=2)
  {
    if (m_dataMode==MODE_HISTOGRAM) Cout("iDataFilePET::PrepareDataFile() -> Build histogram events" << endl);
    else if (m_dataMode==MODE_LIST) Cout("iDataFilePET::PrepareDataFile() -> Build listmode events" << endl);
    else if (m_dataMode==MODE_NORMALIZATION) Cout("iDataFilePET::PrepareDataFile() -> Build normalization events" << endl);
  }

  // ==============================================================================
  // Allocate event buffers (one for each thread)
  // ==============================================================================
  if (m_verbose>=3) Cout("  --> Allocating an event buffer for each thread" << endl);
  // Instanciation of the event buffer according to the data type
  m2p_BufferEvent = new vEvent*[mp_ID->GetNbThreadsForProjection()];
  
  // Allocate the events per each thread
  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    
    // For MODE_LIST events
    if (m_dataMode == MODE_LIST) 
    {
      m2p_BufferEvent[th] = new iEventListPET();
    }
    // For MODE_HISTOGRAM events
    else if (m_dataMode == MODE_HISTOGRAM)
    {
      m2p_BufferEvent[th] = new iEventHistoPET();
      // If we ignore the TOF information, then the event buffer will have only one TOF bin;
      // the TOF contributions will be sum up when reading the event.
      if (m_ignoreTOFFlag) ((iEventHistoPET*)m2p_BufferEvent[th])->SetEventNbTOFBins(1);
      else ((iEventHistoPET*)m2p_BufferEvent[th])->SetEventNbTOFBins(m_nbBinsTOF);
    }
    // For MODE_NORMALIZATION events
    else if (m_dataMode == MODE_NORMALIZATION)
    {
      m2p_BufferEvent[th] = new iEventNorm();
    }

    // Set the maximum number of lines per event
    m2p_BufferEvent[th]->SetNbLines(m_maxNumberOfLinesPerEvent);
    
    // Allocate crystal IDs
    if (m2p_BufferEvent[th]->AllocateID())
    {
      Cerr("*****iDataFilePET::PrepareDataFile() -> Error while trying to allocate memory for the Event object for thread " << th << " !" << endl;);
      return 1;
    }
  }


  // ==============================================================================
  // Deal with specific corrections
  // ==============================================================================

  // In case of attenuation correction flag, see if we ignore this correction
  if (m_atnCorrectionFlag)
  {
    // Affect the ignored flag from the ignored corrections list processed by the oImageDimensionsAndQuantification
    m_ignoreAttnCorrectionFlag = mp_ID->GetIgnoreAttnCorrectionFlag();
    // Verbose
    if (m_verbose>=2)
    {
      if (m_ignoreAttnCorrectionFlag) Cout("  --> Ignore attenuation correction" << endl);
      else Cout("  --> Correct for attenuation" << endl);
    }
  }
  // In case of normalization correction flag, see if we ignore this correction
  if (m_normCorrectionFlag)
  {
    // Affect the ignored flag from the ignored corrections list processed by the oImageDimensionsAndQuantification
    m_ignoreNormCorrectionFlag = mp_ID->GetIgnoreNormCorrectionFlag();
    // Verbose
    if (m_verbose>=2)
    {
      if (m_ignoreNormCorrectionFlag) Cout("  --> Ignore normalization correction" << endl);
      else Cout("  --> Correct for normalization" << endl);
    }
  }
  // In case of scatter correction flag, see if we ignore this correction
  if (m_scatCorrectionFlag)
  {
    // Affect the ignored flag from the ignored corrections list processed by the oImageDimensionsAndQuantification
    m_ignoreScatCorrectionFlag = mp_ID->GetIgnoreScatCorrectionFlag();
    // Verbose
    if (m_verbose>=2)
    {
      if (m_ignoreScatCorrectionFlag) Cout("  --> Ignore scatter correction" << endl);
      else Cout("  --> Correct for scatter events" << endl);
    }
  }
  // In case of random correction flag, see if we ignore this correction
  if (m_randCorrectionFlag)
  {
    // Affect the ignored flag from the ignored corrections list processed by the oImageDimensionsAndQuantification
    m_ignoreRandCorrectionFlag = mp_ID->GetIgnoreRandCorrectionFlag();
    // Verbose
    if (m_verbose>=2)
    {
      if (m_ignoreRandCorrectionFlag) Cout("  --> Ignore random correction" << endl);
      else Cout("  --> Correct for random events" << endl);
    }
  }

  // Normal end
  return 0;
}
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetEventFromBuffer(char* ap_buffer, int a_th)
  \param ap_buffer : address pointing to the event to recover
  \param a_th : index of the thread from which the function was called
  \brief Read an event from the position pointed by 'ap_buffer', parse the generic or modality-specific information, and store them in the (multithreaded) 'm2p_BufferEvent' object
  \return the thread-specific  'm2p_BufferEvent' object containing the modality-specific information for the event
*/
vEvent* iDataFilePET::GetEventFromBuffer(char* ap_buffer, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFilePET::GetEventFromBuffer() ..." << endl);
  #endif

  // Work on a copy of the input pointer
  char* file_position = ap_buffer;

  // For MODE_LIST PET data
  if (m_dataMode == MODE_LIST)
  {
    // Cast the event pointer
    iEventListPET* event = (dynamic_cast<iEventListPET*>(m2p_BufferEvent[a_th]));
    // Mandatory time field: [uint32_t (time)]
    event->SetTimeInMs(*reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // Optional attenuation correction field: [FLTNBDATA (attnCorrectionFactor)]
    if (m_atnCorrectionFlag)
    {
      if (!m_ignoreAttnCorrectionFlag) event->SetAttenuationCorrectionFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }
    // Optional kind: [uint8_t kind]
    if (m_eventKindFlag)
    {
      event->SetKind(*reinterpret_cast<uint8_t*>(file_position));
      file_position += sizeof(uint8_t);
    }
    // Optional scatter correction field: [FLTNBDATA (scatter)]
    if (m_scatCorrectionFlag)
    {
      if (!m_ignoreScatCorrectionFlag) event->SetScatterRate(0, *reinterpret_cast<FLTNBDATA*>(file_position)); 
      file_position += sizeof(FLTNBDATA);
    }
    // Optional random correction field: [FLTNBDATA (random)]
    if (m_randCorrectionFlag)
    {
      if (!m_ignoreRandCorrectionFlag) event->SetRandomRate(*reinterpret_cast<FLTNBDATA*>(file_position)); 
      file_position += sizeof(FLTNBDATA);
    }
    // Optional normalization correction field: [FLTNBDATA (norm)]
    if (m_normCorrectionFlag)
    {
      if (!m_ignoreNormCorrectionFlag) event->SetNormalizationFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }
    // Optional TOF correction field: [uint32_t (TOF)]
    if (m_TOFInfoFlag)
    {
      if (!m_ignoreTOFFlag) event->SetTOFMeasurement(*reinterpret_cast<uint32_t*>(file_position)); 
      file_position += sizeof(uint32_t);
    }
    // Optional POI correction fields: [FLTNBDATA (POI1[3])] [FLTNBDATA (POI2[3])]
    for (int i=0; i<3; i++)
    {
      if (mp_POIDirectionFlag[i])
      {
        if (!m_ignorePOIFlag) event->SetPOI1(i,*reinterpret_cast<FLTNBDATA*>(file_position)); 
        file_position += sizeof(FLTNBDATA);
        if (!m_ignorePOIFlag) event->SetPOI2(i,*reinterpret_cast<FLTNBDATA*>(file_position)); 
        file_position += sizeof(FLTNBDATA);
      }
    }
    // Mandatory nb contributing LORs if m_maxNumberOfLinesPerEvent>1: [uint16_t nbLines]
    if (m_maxNumberOfLinesPerEvent>1)
    {
      event->SetNbLines(*reinterpret_cast<uint16_t*>(file_position));
      file_position += sizeof(uint16_t);
    }
    // Mandatory crystal IDs: [uint32_t (ID1)] [uint32_t (ID2)]
    for (int i=0 ; i<event->GetNbLines() ; i++)
    {
      event->SetID1(i, *reinterpret_cast<uint32_t*>(file_position)); 
      file_position += sizeof(uint32_t);
      event->SetID2(i, *reinterpret_cast<uint32_t*>(file_position)); 
      file_position += sizeof(uint32_t);
    }
  }

  // For MODE_HISTOGRAM PET data
  else if (m_dataMode == MODE_HISTOGRAM) 
  {
    // Cast the event pointer
    iEventHistoPET* event = (dynamic_cast<iEventHistoPET*>(m2p_BufferEvent[a_th]));
    // Mandatory time field: [uint32_t (time)]
    event->SetTimeInMs(*reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // Optional attenuation correction field: [FLTNBDATA (attnCorrectionFactor)]
    if (m_atnCorrectionFlag)
    {
      if (!m_ignoreAttnCorrectionFlag) event->SetAttenuationCorrectionFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }

    // Optional random correction field: [FLTNBDATA (random)]
    if (m_randCorrectionFlag)
    {
      if (!m_ignoreRandCorrectionFlag) event->SetRandomRate(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }

    // Optional normalization correction field: [FLTNBDATA (norm)]
    if (m_normCorrectionFlag)
    {
      if (!m_ignoreNormCorrectionFlag) event->SetNormalizationFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }
    // If we ignore the TOF correction, then we have to sum up the bins' contributions
    if (m_ignoreTOFFlag)
    {
      // Compute the total of event values and scatter rates
      FLTNBDATA total_event_value = 0.;
      FLTNBDATA total_scatter_rate = 0.;
      for (int tb=0 ; tb<m_nbBinsTOF ; tb++)
      {
        // Mandatory bin value: [FLTNBDATA bin value]
        total_event_value += *reinterpret_cast<FLTNBDATA*>(file_position);
        file_position += sizeof(FLTNBDATA);
        // Optional scatter correction field: [FLTNBDATA (scatter)]
        if (m_scatCorrectionFlag)
        {
          total_scatter_rate += *reinterpret_cast<FLTNBDATA*>(file_position);
          file_position += sizeof(FLTNBDATA);
        }
      }
      // Affect the total to the event
      event->SetEventValue(0,total_event_value);
      if (!m_ignoreScatCorrectionFlag) event->SetScatterRate(0,total_scatter_rate);
    }
    // Otherwise we set all different TOF bin contributions
    else
    {
      for (int tb=0 ; tb<m_nbBinsTOF ; tb++)
      {
        // Mandatory bin value: [FLTNBDATA bin value]
        event->SetEventValue(tb, *reinterpret_cast<FLTNBDATA*>(file_position));
        file_position += sizeof(FLTNBDATA);
        // Optional scatter correction field: [FLTNBDATA (scatter)]
        if (m_scatCorrectionFlag)
        {
          if (!m_ignoreScatCorrectionFlag) event->SetScatterRate(tb, *reinterpret_cast<FLTNBDATA*>(file_position));
          file_position += sizeof(FLTNBDATA);
        }
      }
    }
    // Mandatory nb contributing LORs if m_maxNumberOfLinesPerEvent>1: [uint16_t nbLines]
    if (m_maxNumberOfLinesPerEvent>1)
    {
      event->SetNbLines(*reinterpret_cast<uint16_t*>(file_position));
      file_position += sizeof(uint16_t);
    }
    // Mandatory crystal IDs: [uint32_t (c1)] [uint32_t (c2)]
    for (int i=0 ; i<event->GetNbLines() ; i++)
    {
      event->SetID1(i, *reinterpret_cast<uint32_t*>(file_position));
      file_position += sizeof(uint32_t);
      event->SetID2(i, *reinterpret_cast<uint32_t*>(file_position));
      file_position += sizeof(uint32_t); 
    }
  }

  // For MODE_NORMALIZATION PET data
  else if (m_dataMode == MODE_NORMALIZATION)
  {
    // Cast the event pointer
    iEventNorm* event = (dynamic_cast<iEventNorm*>(m2p_BufferEvent[a_th]));
    // Optional attenuation correction field: [FLTNBDATA (norm)]
    if (m_atnCorrectionFlag)
    {
      if (!m_ignoreAttnCorrectionFlag) event->SetAttenuationCorrectionFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }
    // Optional normalization correction field: [FLTNBDATA (norm)]
    if (m_normCorrectionFlag)
    {
      if (!m_ignoreNormCorrectionFlag) event->SetNormalizationFactor(*reinterpret_cast<FLTNBDATA*>(file_position));
      file_position += sizeof(FLTNBDATA);
    }
    // Mandatory nb contributing LORs if m_maxNumberOfLinesPerEvent>1: [uint16_t nbLines]
    if (m_maxNumberOfLinesPerEvent>1)
    {
      event->SetNbLines(*reinterpret_cast<uint16_t*>(file_position));
      file_position += sizeof(uint16_t);
    }
    // Mandatory crystal IDs: [uint32_t (c1)] [uint32_t (c2)]
    for (int i=0 ; i<event->GetNbLines() ; i++)
    {
      event->SetID1(i, *reinterpret_cast<uint32_t*>(file_position));
      file_position += sizeof(uint32_t);
      event->SetID2(i, *reinterpret_cast<uint32_t*>(file_position));
      file_position += sizeof(uint32_t); 
    }
  }

  // Return the updated event
  return m2p_BufferEvent[a_th];
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckSpecificParameters()
  \brief Check parameters specific to PET data
  \todo Adapt m_maxRingDiff initialization for scanner with several layers of crystals (not necessary the same nb of rings)
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::CheckSpecificParameters()
{
  if(m_verbose>=3) Cout("iDataFilePET::CheckSpecificParameters() ..." << endl);

  sScannerManager* p_scannermanager;
  p_scannermanager = sScannerManager::GetInstance(); 
  
  // Error if m_dataType != PET
  if (m_dataType != TYPE_PET)
  {
    Cerr("***** iDataFilePET::CheckSpecificParameters() -> Error, Data type should be PET (=0) !'" << endl);
    return 1;
  }
  
  // Check if Maximum ring difference has been initialized, use the default value from the scanner otherwise.
  if (m_maxRingDiff < 0)
  {
    if (m_verbose>3) Cout("iDataFilePET::CheckSpecificParameters() -> Initialize maxRingDiff with the scanner default value" << endl);
    m_maxRingDiff = p_scannermanager->GetScannerLayerNbRings(0); //TODO maybe throw error here if returned value is incorrect 
                                                                 //TODO  how to deal with scanner with different layers (possibly different number of rings as well)
  }
  
  // Some checks for POI use
  if (m_POIInfoFlag)
  {
    // Not compatible with histogram data
    if (m_dataMode==MODE_HISTOGRAM)
    {
      Cerr("***** iDataFilePET::CheckSpecificParameters() -> POI correction flag is enabled while the data are histogrammed, no sense !" << endl);
      return 1;
    }
    // For each direction (radial, tangential and axial), look at the resolution. If not negative, the info should be the datafile
    for (int i=0; i<3; i++) if (mp_POIResolution[i]>=0.) mp_POIDirectionFlag[i] = true;
  }
  
  // Some checks for TOF use
  if (m_TOFInfoFlag)
  {
    // Check if the resolution was provided
    if (m_resolutionTOF<0.)
    {
      Cerr("***** iDataFilePET::CheckSpecificParameters() -> TOF correction is on while there is no TOF resolution specified in the datafile header !" << endl);
      return 1;
    }
    // For histogram data
    if (m_dataMode==MODE_HISTOGRAM)
    {
      // Check if the number of bins has been provided
      if (m_nbBinsTOF<=1)
      {
        Cerr("***** iDataFilePET::CheckSpecificParameters() -> TOF correction is on while there is only one TOF bin (specified in the datafile header) !" << endl);
        return 1;
      }
      // Check if the bin size has been provided
      if (m_binSizeTOF<=0.)
      {
        Cerr("***** iDataFilePET::CheckSpecificParameters() -> TOF correction is on while there is no bin size specified in the datafile header !" << endl);
        return 1;
      }
    }
  }
  else
  {
    // Set the number of TOF bins to 1
    m_nbBinsTOF = 1;
  }
  
  // End
  return 0;
}
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckFileConsistency()
  \brief  This function is implemented in child classes
          Check if file size is consistent.
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::CheckFileSizeConsistency()
{
  if (m_verbose >=3) Cout("iDataFilePET::CheckFileSizeConsistency() ..." << endl);
    
  // CHECK IF DATAFILE SIZE CORRESPONDS TO THE USER DATA INITIALIZATION 
  m2p_dataFile[0]->seekg(0, ios::end);
  int64_t sizeInBytes = m2p_dataFile[0]->tellg();

  if (m_totalNbEvents*m_sizeEvent != sizeInBytes)
  {
    Cerr("---------------------------------------------------------------------------------------------------------------------------------------------" << endl);
    Cerr("***** iDataFilePET::CheckFileSizeConsistency() -> Error : Datafile size is not consistent with the information provided by the user/datafile!" << endl);
    Cerr("***** iDataFilePET::CheckFileSizeConsistency() -> Expected size : "<< m_totalNbEvents*m_sizeEvent << endl);
    Cerr("***** iDataFilePET::CheckFileSizeConsistency() -> Actual size : "<< sizeInBytes << endl << endl);

    // TODO : Perhaps always write the following depending on verbose (feedback for the user)
    Cerr("***** iDataFilePET::CheckFileSizeConsistency() -> ADDITIONAL INFORMATION ABOUT THE DATAFILE INITIALIZATION : " << endl);

    if(m_maxNumberOfLinesPerEvent > 1)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Compression enabled. Each event should contain " << m_maxNumberOfLinesPerEvent << "LORs, or an equivalent number of LOR + garbage LORs" << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No compression in the data (1 LOR by event in the data) " << endl);
      
    if(m_eventKindFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Coincidence kind (trues, scatter, random, ...) term is enabled " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No informations about the kind of coincidences in the data " << endl);
      
    if(m_normCorrectionFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Normalization correction term is enabled " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No normalization correction terms in the data " << endl);
      
    if(m_atnCorrectionFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Attenuation correction term is enabled " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No attenuation correction terms in the data " << endl);
      
    if(m_scatCorrectionFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Scatter correction term is enabled " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No scatter correction terms in the data " << endl);

    if(m_randCorrectionFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Random correction term is enabled " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No random correction terms in the data " << endl);

    if(m_TOFInfoFlag)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> TOF correction is enabled. Resolution is : " << m_resolutionTOF << " ps."<< endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No TOF correction in the data " << endl);

    Cerr("      iDataFilePET::CheckFileSizeConsistency() -> Calibration factor value is : " << m_calibrationFactor << endl);

    if(mp_POIResolution[0]<0)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No POI enabled on the X axis " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> POI resolution on the X axis is : " << mp_POIResolution[0] << endl);

    if(mp_POIResolution[1]<0)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No POI enabled on the Y axis " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> POI resolution on the Y axis is : " << mp_POIResolution[1] << endl);

    if(mp_POIResolution[2]<0)
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> No POI enabled on the Z axis " << endl);
    else
      Cerr("      iDataFilePET::CheckFileSizeConsistency() -> POI resolution on the Z axis is : " << mp_POIResolution[2] << endl);
    Cerr("---------------------------------------------------------------------------------------------------------------------------------------------" << endl);
    
    return 1;
  }
  
  // Depending on the verbosity, output all the datafile initialization information as feedback for the user
  if (m_verbose>=4)
  {
    Cout("----------------------------------------------------------------Datafile initialization -------------------------------------------------------------" << endl);

    if(m_maxNumberOfLinesPerEvent > 1)
      Cout("  --> Compression enabled. Each event should contain " << m_maxNumberOfLinesPerEvent << "LORs, or an equivalent number of LOR + garbage LORs" << endl);
    else
      Cout("  --> No compression in the data (1 LOR by event in the data) " << endl);

    if(m_eventKindFlag)
      Cout("  --> Coincidence kind (trues, scatter, random, ...) term is enabled " << endl);
    else
      Cout("  --> No informations about the kind of coincidences in the data " << endl);

    if(m_normCorrectionFlag)
      Cout("  --> Normalization correction term is enabled " << endl);
    else
      Cout("  --> No normalization correction terms in the data " << endl);

    if(m_atnCorrectionFlag)
      Cout("  --> Attenuation correction term is enabled " << endl);
    else
      Cout("  --> No attenuation correction terms in the data " << endl);

    if(m_scatCorrectionFlag)
      Cout("  --> Scatter correction term is enabled " << endl);
    else
      Cout("  --> No scatter correction terms in the data " << endl);

    if(m_randCorrectionFlag)
      Cout("  --> Random correction term is enabled " << endl);
    else
      Cout("  --> No random correction terms in the data " << endl);

    if(m_TOFInfoFlag)
      Cout("  --> TOF correction is enabled. Resolution is : " << m_resolutionTOF << " ps."<< endl);
    else
      Cout("  --> No TOF correction in the data " << endl);

    Cout("  --> Calibration factor value is : " << m_calibrationFactor << endl);

    if(mp_POIResolution[0]<0)
      Cout("  --> No POI enabled on the X axis " << endl);
    else
      Cout("  --> POI resolution on the X axis is : " << mp_POIResolution[0] << endl);

    if(mp_POIResolution[1]<0)
      Cout("  --> No POI enabled on the Y axis " << endl);
    else
      Cout("  --> POI resolution on the Y axis is : " << mp_POIResolution[1] << endl);

    if(mp_POIResolution[2]<0)
      Cout("  --> No POI enabled on the Z axis " << endl);
    else
      Cout("  --> POI resolution on the Z axis is : " << mp_POIResolution[2] << endl);
    Cout("---------------------------------------------------------------------------------------------------------------------------------------------" << endl << endl);
  }
  
  // End
  return 0;
}




//PROJECTION SCRIPT FUNCTIONS
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================


/*
  \fn PROJ_InitFile()
  \brief Initialize the fstream objets for output writing as well as some other variables specific to the Projection script (Event-based correction flags, Estimated size of data file)
  \todo Adapt m_maxRingDiff initialization for scanner with several layers of crystals (not necessary the same nb of rings)
  \todo warn the user a datafile with the same name will be erased (eventually add an option to disable the warning)
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::PROJ_InitFile()
{
  if(m_verbose>=3) Cout("iDataFilePET::PROJ_InitFile() ..."<< endl); 
    
  m_startTimeInSec = 0.; //Std initialization for projection
  m_durationInSec = 1.; //Std initialization for projection
  m_totalNbEvents = 0; //Std initialization for projection
  m_nbBinsTOF = 1; // Initialization of this variable required for PROJ_Write functions
  m_isotope = "unknown";

  // Instanciate a fstream datafile for each thread
  m2p_dataFile = new fstream*[mp_ID->GetNbThreadsForProjection()];

  // Set name of the projection data header
  sOutputManager* p_OutMgr;
  p_OutMgr = sOutputManager::GetInstance();
  string path_name = p_OutMgr->GetPathName();
  string img_name = p_OutMgr->GetBaseName();
  m_headerFileName = path_name.append(img_name).append("_CstrProj").append(".Cdh");

  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    m_dataFileName = m_headerFileName.substr(0, m_headerFileName.find_last_of(".")).append(".Cdf"); // Generate datafile name from header file
    
    // We'll write projeted data in several files (corresponding to the number of thread) to be concatenated at the end of the projection process.   
    stringstream ss;
    ss << th;
    string datafile_name = m_dataFileName;
    datafile_name.append("_").append(ss.str());
    
    m2p_dataFile[th] = new fstream( datafile_name.c_str(), ios::binary | ios::out | ios::trunc);
  }


  // Check there is no existing datafile with such name. Remove it otherwise
  ifstream fcheck(m_dataFileName.c_str());
  if(fcheck.good())
  {
    fcheck.close();
    #ifdef _WIN32
    string dos_instruction = "del " + m_dataFileName;
    system(dos_instruction.c_str());
    #else
    remove(m_dataFileName.c_str());
    #endif
  }
  
  if(m_verbose>=3) Cout("iDataFilePET::PROJ_InitFile()-> output datafile name :" << m_dataFileName << endl); 
    
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteEvent(vEvent* ap_Event, FLTNB fp, int a_th, uint32_t a_time)
  \param ap_Event : event containing the data to write
  \param a_th
  \brief  Write event according to the chosen type of data
  \todo   Depending of the RAM load FLAG, either write the data in *nbThreads* different files which will be concatenated at the end (current implementation), 
          or write data in buffers, to be flushed at the end of projection loop.
  \todo   Projection for list-mode
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::PROJ_WriteEvent(vEvent* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFilePET::PROJ_WriteEvent() ..."<< endl); 
  #endif
  
  if(m_dataMode == MODE_LIST) 
  {
    if(PROJ_WriteListEvent((iEventListPET*)ap_Event, a_th))
    {
      Cerr("*****iDataFilePET::PROJ_WriteEvent() -> Error while trying to write projection datafile (list-mode)" << endl;);
      return 1;
    }
  }

  if(m_dataMode == MODE_HISTOGRAM)
  {
    if(PROJ_WriteHistoEvent((iEventHistoPET*)ap_Event, a_th))
    {
      Cerr("*****iDataFilePET::PROJ_WriteEvent() -> Error while trying to write projection datafile (histogram)" << endl;);
      return 1;
    }
  }

  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteHistoEvent()
  \param ap_Event : event containing the data to write
  \param a_th : index of the thread from which the function was called
  \brief  Write a PET histogram event
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::PROJ_WriteHistoEvent(iEventHistoPET* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFilePET::PROJ_WriteHistoEvent() ..."<< endl);
  #endif

  // Write sequentially each field of the event according to the type of the event.
  m2p_dataFile[a_th]->clear();
  
  uint32_t time = ap_Event->GetTimeInMs();
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&time), sizeof(uint32_t));
  
  if(m_atnCorrectionFlag)
  {
    FLTNB atn_corr_factor = ap_Event->GetAtnCorrFactor(); 
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&atn_corr_factor), sizeof(FLTNB));
  }
  
  if(m_randCorrectionFlag) 
  {
    FLTNB rdm_rate = ap_Event->GetEventRdmRate();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&rdm_rate), sizeof(FLTNB));
  }
  
  if(m_normCorrectionFlag)   
  {
    FLTNB norm_corr_factor = ap_Event->GetNormFactor();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&norm_corr_factor), sizeof(FLTNB));
  }
  
  for(int b=0 ; b<m_nbBinsTOF ; b++)
  {
    FLTNB event_value = ap_Event->GetEventValue(b);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&event_value), sizeof(FLTNB));
    if(m_scatCorrectionFlag) 
    {
      FLTNB scat_rate = ap_Event->GetEventScatRate(b);
      m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&scat_rate), sizeof(FLTNB));
    }
  }
  
  uint16_t nb_lines = ap_Event->GetNbLines();
  // Write the number of lines only if the m_maxNumberOfLinesPerEvent is above 1
  if (m_maxNumberOfLinesPerEvent>1) m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&nb_lines), sizeof(uint16_t));

  for(int i=0 ; i<nb_lines ; i++)
  {
    uint32_t id1 = ap_Event->GetID1(i);
    uint32_t id2 = ap_Event->GetID2(i);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id1), sizeof(uint32_t));
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id2), sizeof(uint32_t));
  }

  // 0-filling if needed (nb lines inferior to the max number for PET data with compression)
  if(nb_lines<m_maxNumberOfLinesPerEvent)
    for(int i=0 ; i<m_maxNumberOfLinesPerEvent-nb_lines ; i++)
    {
      uint32_t gbg = 0;
      m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&gbg), sizeof(uint32_t));
      m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&gbg), sizeof(uint32_t));
    }

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteListEvent()
  \param ap_Event : event containing the data to write
  \param a_th : index of the thread from which the function was called
  \brief  Write a PET list-mode event
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::PROJ_WriteListEvent(iEventListPET* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFilePET::PROJ_WriteListEvent() ..."<< endl);
  #endif

  // Write sequentially each field of the event according to the type of the event.
  m2p_dataFile[a_th]->clear();
  
  uint32_t time = ap_Event->GetTimeInMs();
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&time), sizeof(uint32_t));
  
  
  if(m_atnCorrectionFlag)
  {
    FLTNB atn_corr_factor = ap_Event->GetAtnCorrFactor(); 
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&atn_corr_factor), sizeof(FLTNB));
  }
  
  if(m_eventKindFlag) 
  {
    uint8_t event_kind = ap_Event->GetKind();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&event_kind), sizeof(uint8_t));
  }
  
  if(m_scatCorrectionFlag)   
  {
    FLTNB scat_rate = ap_Event->GetEventScatRate(0);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&scat_rate), sizeof(FLTNB));
  }
  
  if(m_randCorrectionFlag) 
  {
    FLTNB rdm_rate = ap_Event->GetEventRdmRate();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&rdm_rate), sizeof(FLTNB));
  }
  
  if(m_normCorrectionFlag)   
  {
    FLTNB norm_corr_factor = ap_Event->GetNormFactor();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&norm_corr_factor), sizeof(FLTNB));
  }
    
  if(m_TOFInfoFlag)
  {
    uint32_t TOF_measurement = ap_Event->GetTOFMeasurement();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&TOF_measurement), sizeof(uint32_t));
  }
  
  
  if(mp_POIResolution[0]>0)
  {
    FLTNB POI_1 = ap_Event->GetPOI1(0);
    FLTNB POI_2 = ap_Event->GetPOI2(0);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_1), sizeof(FLTNB));
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_2), sizeof(FLTNB));
  }
  if(mp_POIResolution[1]>0)
  {
    FLTNB POI_1 = ap_Event->GetPOI1(1);
    FLTNB POI_2 = ap_Event->GetPOI2(1);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_1), sizeof(FLTNB));
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_2), sizeof(FLTNB));
  }
  if(mp_POIResolution[2]>0)
  {
    FLTNB POI_1 = ap_Event->GetPOI1(2);
    FLTNB POI_2 = ap_Event->GetPOI2(2);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_1), sizeof(FLTNB));
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI_2), sizeof(FLTNB));
  }
  
  uint16_t nb_lines = ap_Event->GetNbLines();
  // Write the number of lines only if the m_maxNumberOfLinesPerEvent is above 1
  if (m_maxNumberOfLinesPerEvent>1) m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&nb_lines), sizeof(uint16_t));

  for(int i=0 ; i<nb_lines ; i++)
  {
    uint32_t id1 = ap_Event->GetID1(i);
    uint32_t id2 = ap_Event->GetID2(i);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id1), sizeof(uint32_t));
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id2), sizeof(uint32_t));
  }
  
  // 0-filling if needed (nb lines inferior to the max number for PET data with compression)
  if(nb_lines<m_maxNumberOfLinesPerEvent)
    for(int i=0 ; i<m_maxNumberOfLinesPerEvent-nb_lines ; i++)
    {
      uint32_t gbg = 0;
      m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&gbg), sizeof(uint32_t));
      m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&gbg), sizeof(uint32_t));
    }

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteHeader()
  \brief  Generate a header file according to the projection and data output informations.
          Used by Projection algorithm.
  \todo   Use Interfile format
  \return 0 if success, and positive value otherwise.
*/
int iDataFilePET::PROJ_WriteHeader()
{
  if(m_verbose>=3) Cout("PROJ_WriteHeader() ..."<< endl);
  
  fstream headerFile;
  stringstream ss;

  headerFile.open(m_headerFileName.c_str(), ios::out);

  // ----- MANDATORY FLAGS ----- // 
  headerFile << "Data filename:" << "    " << GetFileFromPath(m_dataFileName).c_str() << endl; 
  headerFile << "Number of events:" << "    " << m_totalNbEvents << endl; 

  //- Flag: list-mode or histogram mode
  headerFile << "Data mode:" << "    " << m_dataMode << endl; 
  // Flag: PET, SPECT or TRANSMISSION
  headerFile << "Data type: " << "    " << "PET" << endl;

  headerFile << "Start time (s):" << "    " << m_startTimeInSec << endl; 
  
  headerFile << "Duration (s):" << "    " << m_durationInSec << endl; 

  sScannerManager* p_scannermanager; 
  p_scannermanager = sScannerManager::GetInstance(); 

  headerFile << "Scanner name:" << "    " << p_scannermanager->GetScannerName() << endl; 
  
  headerFile << "Maximum ring difference:" << "    " << m_maxRingDiff << endl; 
  
  // Only write this field if different than default value
  if(m_maxNumberOfLinesPerEvent >1)
    headerFile << "Maximum number of lines per event:" << "    " << m_maxNumberOfLinesPerEvent << endl; 
    
    
  
  // ----- OPTIONNAL FLAGS ----- // 
  headerFile << "Calibration factor:" << "    " << m_calibrationFactor << endl;

  if (m_isotope != "unknown") headerFile << "Isotope:" << "    " << m_isotope << endl; 
    
    

  // ----- TOF info ----- // 
  if(m_TOFInfoFlag)
  {
    if(m_dataMode == MODE_LIST) 
    {
      if(m_resolutionTOF <= 0.) // Check if TOF resolution (list-mode) initialized, error otherwise
      { 
        Cerr("***** iDataFilePET::PROJ_WriteHeader -> Error : TOF resolution (mm) for list-mode not initialized, while TOF correction flag is enabled!" << endl);
        return 1;
      }
      else
        headerFile << "TOF resolution:" << "    " << m_resolutionTOF << endl; 
    }
    else if (m_dataMode == MODE_HISTOGRAM) 
    {
      if(m_binSizeTOF > 0.) // Check if TOF bin size (histogram) initialized, error otherwise
      { 
        Cerr("***** iDataFilePET::PROJ_WriteHeader -> Error : TOF bin size for histogram mode not initialized, while TOF correction flag is enabled!" << endl);
        return 1;
      }
      else
      {
        headerFile << "Histo TOF bin:" << "    " << m_nbBinsTOF << endl; 
        headerFile << "Histo TOF size:" << "    " << m_binSizeTOF << endl; 
      }
    }
    else
    { 
      Cerr("***** iDataFilePET::PROJ_WriteHeader -> Error : Datafile mode unknown (should be 0 (list-mode) or histogram (1))!" << endl);
      return 1;
    }
  }
  
  
  // ----- POI info ----- // 
  if (mp_POIResolution[0]>0. || mp_POIResolution[1]>0. || mp_POIResolution[2]>0.)
  {
    // Error if histogram data are used
    if (m_dataMode==MODE_HISTOGRAM)
    {
      Cerr("***** iDataFilePET::PROJ_WriteHeader -> Request for writing POI resolution in histogram mode (this field is specific to list-mode) !" << endl);
      return 1;
    }
    else
      headerFile << "DOI capability:" << "    " << mp_POIResolution[0] << "," << mp_POIResolution[1] << "," << mp_POIResolution[2] << endl;
  }


  // ----- Event kind ----- // 
  if (m_eventKindFlag)
  {
    // Error if histogram data are used
    if(m_dataMode==MODE_HISTOGRAM)
    {
      Cerr("***** iDataFilePET::PROJ_WriteHeader -> Request for writing event type in histogram mode (this field is specific to list-mode) !" << endl);
      return 1;
    }
    else 
      headerFile << "Coincidence kind informations flag:" << "    " << m_eventKindFlag << endl;
  }
 
 
   // ----- Correction flags ----- // 
   
  if(m_atnCorrectionFlag)
    headerFile << "Attenuation correction flag:" << "    " << m_atnCorrectionFlag << endl;  //TODO Should be computed during projection
    
  if(m_normCorrectionFlag)
    headerFile << "Normalization correction flag:" << "    " << m_normCorrectionFlag << endl; 

  if(m_scatCorrectionFlag)
    headerFile << "Scatter correction flag:" << "    " << m_scatCorrectionFlag << endl; 
    
  if(m_randCorrectionFlag)
    headerFile << "Random correction flag:" << "    " << m_randCorrectionFlag << endl; 
  
  headerFile.close();

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_GetScannerSpecificParameters()
  \brief Get PET specific parameters for projections from the scanner object, through the scannerManager.
  \return 0 if success, positive value otherwise
*/
int iDataFilePET::PROJ_GetScannerSpecificParameters()
{
  if(m_verbose>=3) Cout("PROJ_GetScannerSpecificParameters() ..."<< endl);
  sScannerManager* p_scannermanager;
  p_scannermanager = sScannerManager::GetInstance(); 
   
  // Get pointers to SPECT specific parameters in scanner
  if(p_scannermanager->PROJ_GetPETSpecificParameters(&m_maxRingDiff) )
  {
    Cerr("***** iDataFilePET::PROJ_GetScannerSpecificParameters() -> An error occurred while trying to get PET geometric parameters from the scanner object !" << endl);
    return 1;
  }

  return 0;
}