iDataFileSPECT.cc

Go to the documentation of this file.

#include "iDataFileSPECT.hh"

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

iDataFileSPECT::iDataFileSPECT() : vDataFile() 
{
  // Set all members to default values
  m_dataType = TYPE_SPECT;
  m_isotope = "unknown";
  mp_nbOfBins[0] = 1;
  mp_nbOfBins[1] = 1;
  m_nbOfProjections = 0;
  mp_angles = NULL;
  m_nbHeads = 1;
  mp_CORtoDetectorDistance = NULL;
  m_eventKindFlag = false;
  m_normCorrectionFlag = false;
  m_ignoreNormCorrectionFlag = false;
  m_scatCorrectionFlag = false;
  m_ignoreScatCorrectionFlag = false;
  m_headRotDirection = GEO_ROT_CW;
}

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

iDataFileSPECT::~iDataFileSPECT() 
{
  if (mp_angles) delete[] mp_angles;
  if (mp_CORtoDetectorDistance) delete[] mp_CORtoDetectorDistance;
}

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

int iDataFileSPECT::ReadSpecificInfoInHeader(bool a_affectQuantificationFlag)
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Verbose
  if (m_verbose>=VERBOSE_DETAIL) Cout("iDataFileSPECT::ReadSpecificInfoInHeader() -> Read information specific to SPECT" << endl);

  // Create pointers dedicated to recover the addresses of the member variables of the scanner object
  FLTNB* p_angles = NULL;
  FLTNB* p_CORtoDetectorDistance = NULL; 
  FLTNB p_pixSizeXY[2]; // not used here
  
  // Get Geometric parameters recovered from the scanner object
  sScannerManager* p_scannermanager;
  p_scannermanager = sScannerManager::GetInstance(); 
  p_scannermanager->GetSPECTSpecificParameters(&m_nbOfProjections, 
                                               &m_nbHeads, 
                                                mp_nbOfBins, 
                                                p_pixSizeXY, 
                                                p_angles, 
                                                p_CORtoDetectorDistance,
                                               &m_headRotDirection);

  // Check m_nbOfProjections first before allocating projection angles and radius using this variable
  if (m_nbOfProjections==0)
  {
    Cerr("***** iDataFileSPECT::ReadSpecificInfoInHeader() -> Number of projections should be strictly positive !" << endl);
    return 1;
  }

  // Allocation and initialization of Projection angles and Center of rotation to SPECT detector surface distance (radius) :
  mp_angles = new FLTNB[m_nbOfProjections];
  mp_CORtoDetectorDistance = new FLTNB[m_nbOfProjections];

  // Recover values
  for (int a=0 ; a<m_nbOfProjections ; a++)
  {
    mp_angles[a] = p_angles[a];
    mp_CORtoDetectorDistance[a] = p_CORtoDetectorDistance[a];
  }

  // Feedback to user
  if (m_verbose==VERBOSE_DETAIL) 
  {
    Cout("  --> Provided projection angles | distance Center of Rotation (COR) to detector: " << endl);
    for (int a=0 ; a<m_nbOfProjections ; a++) Cout("      " << mp_angles[a] << " | "  << mp_CORtoDetectorDistance[a] << endl);
  }

  // Read optional fields in the header, check if errors (issue during data reading/conversion (==1) )
  if (ReadDataASCIIFile(m_headerFileName, "Isotope", &m_isotope, 1, 0) == 1 ||
      ReadDataASCIIFile(m_headerFileName, "Event kind flag", &m_eventKindFlag, 1, 0) == 1 ||
      ReadDataASCIIFile(m_headerFileName, "Normalization correction flag", &m_normCorrectionFlag, 1, 0) == 1 ||
      ReadDataASCIIFile(m_headerFileName, "Scatter correction flag", &m_scatCorrectionFlag, 1, 0) == 1 )
      {
        Cerr("***** iDataFileSPECT::ReadSpecificInfoInHeader() -> Error while reading optional fields in the header data file !" << endl);
        return 1;
      }

  // Give the SPECT isotope to the oImageDimensionsAndQuantification that manages the quantification factors
  if (a_affectQuantificationFlag && mp_ID->SetSPECTIsotope(m_bedIndex, m_isotope))
  {
    Cerr("***** iDataFileSPECT::ReadSpecificInfoInHeader() -> A problem occured while setting the isotope to oImageDimensionsAndQuantification !" << endl);
    return 1;
  }

  // Normal end
  return 0;
}

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

int iDataFileSPECT::ComputeSizeEvent()
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Verbose
  if (m_verbose>=VERBOSE_DETAIL) Cout("iDataFileSPECT::ComputeSizeEvent() -> In bytes" << endl);

  // For MODE_LIST events
  if (m_dataMode == MODE_LIST) 
  {
    // Size of the mandatory element in a list-mode event: Time + 2*crystalID
    m_sizeEvent = sizeof(uint32_t) + 2*sizeof(uint32_t);
    // Optional flags
    if (m_eventKindFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_scatCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_normCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    for (int i=0 ; i<3; i++)
    {
      // POI available for this direction
      if (mp_POIDirectionFlag[i]) m_sizeEvent += sizeof(FLTNBDATA);
    }
  }
  // For MODE_HISTOGRAM events
  if(m_dataMode == MODE_HISTOGRAM)
  {
    // Size of the mandatory element in a histo event: Time + event_value + 2*crystalID
    m_sizeEvent = sizeof(uint32_t) + sizeof(FLTNBDATA) + 2*sizeof(uint32_t);
    // Optional flags
    if (m_scatCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
    if (m_normCorrectionFlag) m_sizeEvent += sizeof(FLTNBDATA);
  }
  // Unknown event type
  else
  {
    Cerr("***** iDataFileSPECT::ComputeSizeEvent() -> Unknown event mode !" << endl);
    return 1;
  }

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

  // Verbose
  if (m_verbose>=VERBOSE_DETAIL) Cout("  --> Event size: " << m_sizeEvent << " bytes" << endl);
  // End
  return 0;
}

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

int iDataFileSPECT::PrepareDataFile()
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Verbose
  if (m_verbose>=VERBOSE_NORMAL)
  {
    if (m_dataMode==MODE_HISTOGRAM) Cout("iDataFileSPECT::PrepareDataFile() -> Build histogram events" << endl);
    else if (m_dataMode==MODE_LIST) Cout("iDataFileSPECT::PrepareDataFile() -> Build listmode events" << endl);
    else if (m_dataMode==MODE_NORMALIZATION) Cout("iDataFileSPECT::PrepareDataFile() -> Build normalization events" << endl);
  }

  // ==============================================================================
  // Allocate event buffers (one for each thread)
  // ==============================================================================

  if (m_verbose>=VERBOSE_DETAIL) Cout("  --> Allocate 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 iEventListSPECT();
    }
    // For MODE_HISTOGRAM events
    if (m_dataMode == MODE_HISTOGRAM)
    {
      m2p_BufferEvent[th] = new iEventHistoSPECT();
    }
    // Allocate crystal/view IDs
    if (m2p_BufferEvent[th]->AllocateID())
    {
      Cerr("*****iDataFileSPECT::PrepareDataFile() -> Error while trying to allocate memory for the Event object!" << endl);
      return 1;
    }
  }

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

  // 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>=VERBOSE_DETAIL)
    {
      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>=VERBOSE_DETAIL)
    {
      if (m_ignoreScatCorrectionFlag) Cout("  --> Ignore scatter correction" << endl);
      else Cout("  --> Correct for scatter events" << endl);
    }
  }

  // Normal end
  return 0;
}

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

vEvent* iDataFileSPECT::GetEventFromBuffer(char* ap_buffer, int a_th)
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_EVENT)

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

  // For MODE_LIST SPECT data
  if (m_dataMode == MODE_LIST)
  {
    // Cast the event pointer
    iEventListSPECT* event = (dynamic_cast<iEventListSPECT*>(m2p_BufferEvent[a_th]));
    // Mandatory time field: [uint32_t (time)]
    event->SetTimeInMs(*reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // 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(*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 POI correction field: [FLTNBDATA (POI1[3])]
    for (int i=0; i<3; i++)
    {
      if (mp_POIDirectionFlag[i])
      {
        if (!m_ignorePOIFlag) event->SetPOI(i,*reinterpret_cast<FLTNBDATA*>(file_position)); 
        file_position += sizeof(FLTNBDATA);
      }
    }
    // Mandatory angular projection ID: [uint32_t (ID1)]
    event->SetID1(0, *reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // Mandatory crystal ID: [uint32_t (ID2)]
    event->SetID2(0, *reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
  }

  // For MODE_HISTOGRAM SPECT DATA
  if (m_dataMode == MODE_HISTOGRAM) 
  {
    // Cast the event pointer
    iEventHistoSPECT* event = (dynamic_cast<iEventHistoSPECT*>(m2p_BufferEvent[a_th]));
    // Mandatory time field: [uint32_t (time)]
    event->SetTimeInMs(*reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // Mandatory bin value: [FLTNBDATA bin value]
    event->SetEventValue(0, *reinterpret_cast<FLTNBDATA*>(file_position));
    file_position += sizeof(FLTNBDATA);
    // Optional scatter correction field: [FLTNBDATA (scatter)]
    if (m_scatCorrectionFlag)
    {
      if (!m_ignoreScatCorrectionFlag) event->SetScatterRate(*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 angular projection ID: [uint32_t (ID1)]
    event->SetID1(0, *reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
    // Mandatory crystal ID: [uint32_t (ID2)]
    event->SetID2(0, *reinterpret_cast<uint32_t*>(file_position)); 
    file_position += sizeof(uint32_t);
  }

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

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

int iDataFileSPECT::CheckSpecificParameters()
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Error if m_dataType != SPECT
  if (m_dataType != TYPE_SPECT)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Data type should be SPECT !'" << endl);
    return 1;
  }
  // Check number of projections
  if (m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Number of projections not initialized (should be >0) !" << endl);
    return 1;
  }
  // Check projection angles
  if (mp_angles == NULL)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Projection angles not initialized !" << endl);
    return 1;
  }
  // End
  return 0;
}

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

int iDataFileSPECT::CheckFileSizeConsistency()
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  
  // Check datafile self-consistency
  m2p_dataFile[0]->seekg(0, ios::end);
  int64_t sizeInBytes = m2p_dataFile[0]->tellg();
  if (m_totalNbEvents*m_sizeEvent != sizeInBytes)
  {
    Cerr("----------------------------------------------------------------------------------------------------------------------------------------" << endl);
    Cerr("***** iDataFileSPECT::CheckFileSizeConsistency() -> Datafile size is not consistent with the information provided by the user/datafile !" << endl);
    Cerr("  --> Expected size: "<< m_totalNbEvents*m_sizeEvent << endl);
    Cerr("  --> Actual size: "<< sizeInBytes << endl << endl);
    Cerr("      ADDITIONAL INFORMATION ABOUT THE DATAFILE INITIALIZATION" << endl);
    if (m_eventKindFlag) Cerr("  --> Coincidence kind term is enabled" << endl);
    else Cerr("  --> No information about the kind of coincidences in the data" << endl);
    if (m_normCorrectionFlag) Cerr("  --> Normalization correction term is enabled" << endl);
    else Cerr("  --> No normalization correction term in the data" << endl);
    if (m_scatCorrectionFlag) Cerr("  --> Scatter correction term is enabled" << endl);
    else  Cerr("  --> No scatter correction term in the data" << endl);
    Cerr("  --> Calibration factor value is: " << m_calibrationFactor << endl);
    Cerr("  --> Number of bins are equal to " << mp_nbOfBins[0] << "," << mp_nbOfBins[1] << " for transaxial,axial axis respectively" << endl);
    if (mp_POIResolution[0]<0.) Cerr("  --> No POI enabled on the radial axis" << endl);
    else Cerr("  --> POI resolution on the radial axis is: " << mp_POIResolution[0] << endl);
    if (mp_POIResolution[1]<0.) Cerr("  --> No POI enabled on the tangential axis" << endl);
    else Cerr("  --> POI resolution on the tangential axis is: " << mp_POIResolution[1] << endl);
    if (mp_POIResolution[2]<0.) Cerr("  --> No POI enabled on the axial axis" << endl);
    else Cerr("  --> POI resolution on the axial 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>=VERBOSE_DETAIL)
  {
    Cout("------------------------------------------------- Datafile initialization -------------------------------------------------" << endl);
    if (m_eventKindFlag) Cout("  --> Coincidence kind term is enabled" << endl);
    else Cout("  --> No information 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_scatCorrectionFlag) Cout("  --> Scatter correction term is enabled" << endl);
    else  Cout("  --> No scatter correction terms in the data" << endl);
    Cout("  --> Calibration factor value is: " << m_calibrationFactor << endl);
    Cout("  --> Number of bins are equal to " << mp_nbOfBins[0] << "," << mp_nbOfBins[1] << " for transaxial,axial axis respectively" << endl);
    if (mp_POIResolution[0]<0.) Cout("  --> No POI enabled on the radial axis" << endl);
    else Cout("  --> POI resolution on the radial axis is: " << mp_POIResolution[0] << endl);
    if (mp_POIResolution[1]<0.) Cout("  --> No POI enabled on the tangential axis" << endl);
    else Cout("  --> POI resolution on the tangential axis is: " << mp_POIResolution[1] << endl);
    if (mp_POIResolution[2]<0.) Cout("  --> No POI enabled on the axial axis" << endl);
    else Cout("  --> POI resolution on the axial axis is: " << mp_POIResolution[2] << endl);
    Cout("---------------------------------------------------------------------------------------------------------------------------" << endl);
  }

  // End
  return 0;
}

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

int iDataFileSPECT::CheckSpecificConsistencyWithAnotherDatafile(vDataFile* ap_Datafile)
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Dynamic cast the vDataFile to a iDataFilePET
  iDataFileSPECT* p_data_file = (dynamic_cast<iDataFileSPECT*>(ap_Datafile));
  // Check isotope
  if (m_isotope!=p_data_file->GetIsotope())
  {
    Cerr("***** iDataFileSPECT::CheckSpecificConsistencyWithAnotherDatafile() -> Isotopes are inconsistent !" << endl);
    return 1;
  }
  // Check event kind flag
  if (m_eventKindFlag!=p_data_file->GetEventKindFlag())
  {
    Cerr("***** iDataFileSPECT::CheckSpecificConsistencyWithAnotherDatafile() -> Event kind flags are inconsistent !" << endl);
    return 1;
  }
  // Check normalization correction flag
  if (m_normCorrectionFlag!=p_data_file->GetNormCorrectionFlag())
  {
    Cerr("***** iDataFileSPECT::CheckSpecificConsistencyWithAnotherDatafile() -> Normalization correction flags are inconsistent !" << endl);
    return 1;
  }
  // Check scatter correction flag
  if (m_scatCorrectionFlag!=p_data_file->GetScatCorrectionFlag())
  {
    Cerr("***** iDataFileSPECT::CheckSpecificConsistencyWithAnotherDatafile() -> Scatter correction flags are inconsistent !" << endl);
    return 1;
  }
  // End
  return 0;
}

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

int iDataFileSPECT::InitAngles(FLTNB* ap_angles)
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Check
  if (m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::InitAngles() -> Number of projection angles not initialized !'" << endl);
    return 1;
  }
  // Allocate
  mp_angles = new FLTNB[m_nbOfProjections];
  // Affect
  for (uint16_t a=0 ; a<m_nbOfProjections ; a++) mp_angles[a] = ap_angles[a];
  // End
  return 0;
}

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

int iDataFileSPECT::InitCorToDetectorDistance(FLTNB* ap_CORtoDetectorDistance)
{
  DEBUG_VERBOSE(m_verbose,VERBOSE_DEBUG_LIGHT)
  // Check
  if (m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::InitAngles() -> Number of projection angles not initialized !'" << endl);
    return 1;
  }
  // Allocate
  mp_CORtoDetectorDistance = new FLTNB[m_nbOfProjections];
  // Affect
  for (uint16_t a=0 ; a<m_nbOfProjections ; a++) mp_CORtoDetectorDistance[a] = ap_CORtoDetectorDistance[a];
  // 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 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 iDataFileSPECT::PROJ_InitFile()
{
  if(m_verbose >=3) Cout("iDataFileSPECT::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_calibrationFactor =  1.;
  m_isotope = "unknown";
  m_scatCorrectionFlag = false;
  m_normCorrectionFlag = false;
  mp_POIResolution[0] = -1.;
  mp_POIResolution[1] = -1.;
  mp_POIResolution[2] = -1.;

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

    // Projeted data will be written 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);
  }

  //remove content from the output data file, in case it already exists
  //todo warn the user a datafile with the same name will be erased (eventually add an option to disable the warning)
  ofstream output_file(m_dataFileName.c_str(), ios::out | ios::trunc);
  output_file.close();
  
  if(m_verbose>=3) Cout("iDataFileSPECT::PROJ_InitFile()-> output datafile name :" << m_dataFileName << endl); 
    
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_WriteEvent(vEvent* ap_Event, int a_th)
  \param ap_Event : event containing the data to write
  \param a_th : index of the thread from which the function was called
  \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 iDataFileSPECT::PROJ_WriteEvent(vEvent* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFileSPECT::PROJ_WriteEvent() ..."<< endl); 
  #endif
  
  if(m_dataMode == MODE_LIST) 
  {
    // TODO should create as many vEvent as (int)fp.
    if(PROJ_WriteListEvent((iEventListSPECT*)ap_Event, a_th))
    {
      Cerr("*****iDataFileSPECT::PROJ_WriteEvent() -> Error while trying to write projection datafile (list-mode)" << endl);
      return 1;
    }
  }

  if(m_dataMode == MODE_HISTOGRAM)
  {
    if(PROJ_WriteHistoEvent((iEventHistoSPECT*)ap_Event, a_th))
    {
      Cerr("*****iDataFileSPECT::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 SPECT histogram event
  \return 0 if success, and positive value otherwise.
*/
int iDataFileSPECT::PROJ_WriteHistoEvent(iEventHistoSPECT* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFileSPECT::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));
  
  FLTNB event_value = ap_Event->GetEventValue(0);
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&event_value), sizeof(FLTNB));

  if(m_scatCorrectionFlag)   
  {
    FLTNB scat_rate = ap_Event->GetEventScatRate();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&scat_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));
  }

  uint32_t id1 = ap_Event->GetID1(0);
  uint32_t id2 = ap_Event->GetID2(0);
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id1), sizeof(uint32_t));
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id2), 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 SPECT list-mode event
  \return 0 if success, and positive value otherwise.
*/
int iDataFileSPECT::PROJ_WriteListEvent(iEventListSPECT* ap_Event, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFileSPECT::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_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();
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&scat_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(mp_POIResolution[0]>0)
  {
    FLTNB POI = ap_Event->GetPOI(0);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI), sizeof(FLTNB));
  }

  if(mp_POIResolution[1]>0)
  {
    FLTNB POI = ap_Event->GetPOI(1);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI), sizeof(FLTNB));
  }

  if(mp_POIResolution[2]>0)
  {
    FLTNB POI = ap_Event->GetPOI(2);
    m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&POI), sizeof(FLTNB));
  }

  uint32_t id1 = ap_Event->GetID1(0);
  uint32_t id2 = ap_Event->GetID2(0);
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id1), sizeof(uint32_t));
  m2p_dataFile[a_th]->write(reinterpret_cast<char*>(&id2), 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.
  \return 0 if success, and positive value otherwise.
*/
int iDataFileSPECT::PROJ_WriteHeader()
{
  // Verbose
  if(m_verbose>=3) Cout("iDataFileSPECT::PROJ_WriteHeader() ..."<< endl); 

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

  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: " << "    " << "SPECT" << 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; 

  
  if(mp_nbOfBins[0] != 0 && mp_nbOfBins[1] != 0)
    headerFile << "Number of bins:" << "    " << mp_nbOfBins[0] << "," << mp_nbOfBins[1] << endl; 
  
  headerFile << "Number of projections:" << "    " << m_nbOfProjections << endl; 

  headerFile << "Projection angles:" << "    "; 
  
  headerFile << mp_angles[0]; 
  for(int a=1 ; a<m_nbOfProjections ; a++)
    headerFile << ", " << mp_angles[a]; 

  headerFile << endl;
  if(m_nbHeads>1)
  {
    headerFile << "Global distance camera surface to COR:" << "    " << mp_CORtoDetectorDistance[0];
    for(int h=1 ; h<m_nbHeads ; h++) headerFile << "," << mp_CORtoDetectorDistance[h] << endl;
    headerFile << endl;
  }
  
  else if(mp_CORtoDetectorDistance[0] > 0)
    headerFile << "Global distance camera surface to COR:" << "    " << mp_CORtoDetectorDistance[0] << endl;
  
  // Optional fields
  headerFile << "Calibration factor:" << "    " << m_calibrationFactor << endl;
  headerFile << "Isotope:" << "    " << m_isotope << endl; 
  headerFile << "DOI capability:" << "    " << mp_POIResolution[0] << "," << mp_POIResolution[1] << "," << mp_POIResolution[2] << endl; 
  headerFile << "Normalization correction flag:" << "    " << m_normCorrectionFlag << endl; 
  headerFile << "Scatter correction flag:" << "    " << m_scatCorrectionFlag << endl; 
  string rot_direction = (m_headRotDirection == GEO_ROT_CCW) ? "CCW" : "CW";
  headerFile << "Head rotation direction:" << "    " << rot_direction << endl; 
  
  headerFile.close();

  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_GetScannerSpecificParameters()
  \brief Get SPECT specific parameters for projections from the scanner object, through the scannerManager.
  \return 0 if success, positive value otherwise
*/
int iDataFileSPECT::PROJ_GetScannerSpecificParameters()
{
  // Verbose
  if(m_verbose>=3) Cout("iDataFileSPECT::PROJ_GetScannerSpecificParameters() ..."<< endl); 
  
  sScannerManager* p_scannermanager;
  p_scannermanager = sScannerManager::GetInstance(); 
   
  // Create pointers dedicated to recover the addresses of the member variables of the scanner object
  FLTNB* p_angles = NULL;
  FLTNB* p_CORtoDetectorDistance = NULL; 
  FLTNB p_pixSizeXY[2] ; // not used here
  
  // Get pointers to SPECT specific parameters in scanner
  if(p_scannermanager->GetSPECTSpecificParameters(&m_nbOfProjections, &m_nbHeads, mp_nbOfBins, p_pixSizeXY, p_angles, p_CORtoDetectorDistance, &m_headRotDirection) )
  {
    Cerr("*****iDataFileSPECT::PROJ_GetScannerSpecificParameters() -> An error occurred while trying to get SPECT geometric parameters from the scanner object !" << endl);
    return 1;
  }

  // In projection, mp_CORtoDetectorDistance is allocated according to the number of heads, and not the number of projections
  // (Projection does not currently allow to set a radius specific to the projection angles)
  mp_CORtoDetectorDistance = new FLTNB[m_nbHeads];

  // Retrieve SPECT distance between the detector and center of rotation
  for(int h=0 ; h<m_nbHeads ; h++)
    mp_CORtoDetectorDistance[h] = p_CORtoDetectorDistance[0];

  // Retrieve projection angles
  mp_angles = new FLTNB[m_nbOfProjections];

  for(int a=0 ; a<m_nbOfProjections ; a++)
    mp_angles[a] = p_angles[a];
      
  return 0;
}