iDataFileSPECT.cc

Go to the documentation of this file.

/*
  Implementation of class iDataFileSPECT

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



#include "iDataFileSPECT.hh"



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iDataFileSPECT constructor. 
         Initialize the member variables to their default values.
*/
iDataFileSPECT::iDataFileSPECT() : vDataFile() 
{
  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;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iDataFileSPECT destructor. 
*/
iDataFileSPECT::~iDataFileSPECT() 
{
  if (mp_angles) delete[] mp_angles;
  if (mp_CORtoDetectorDistance) delete[] mp_CORtoDetectorDistance;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadSpecificInfoInHeader()
  \brief Read through the header file and recover specific SPECT information.
  \return 0 is success, positive value otherwise
*/
int iDataFileSPECT::ReadSpecificInfoInHeader()
{
  // Verbose
  if (m_verbose>=3) Cout("iDataFileSPECT::ReadSpecificInfoInHeader() ..." << 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 mp_nbOfBins[2]
  /*
  if(mp_nbOfBins[0] == 0 || mp_nbOfBins[1] == 0)
  {
    Cerr("***** iDataFileSPECT::ReadSpecificInfoInHeader() -> Error, number of bins should be >0 '" << endl);
    return 1;
  }
*/
  // Check m_nbOfProjections first before allocating projection angles and radius using this variable
  if(m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::ReadSpecificInfoInHeader() -> Error, number of projections should be >0 '" << 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 ==3) 
  {
    Cout("iDataFileSPECT::ReadSpecificInfoInHeader() -> 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);
    Cout(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 field in the header data file '" << endl);
        return 1;
      }

  // Give the SPECT isotope to the oImageDimensionsAndQuantification that manages the quantification factors
  if (m_dataMode!=MODE_NORMALIZATION && mp_ID->SetSPECTIsotope(m_bedIndex, m_isotope))
  {
    Cerr("***** iDataFileSPECT::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
  \return 0 is success, positive value otherwise
*/
int iDataFileSPECT::ComputeSizeEvent()
{
  // Verbose
  if (m_verbose>=3) Cout("iDataFileSPECT::ComputeSizeEvent() ..." << 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);
  }

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

  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 iDataFileSPECT::PrepareDataFile()
{
  // Verbose
  if (m_verbose>=2)
  {
    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>=3) 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>=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);
    }
  }

  // Normal end
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetEventFromBuffer()
  \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* iDataFileSPECT::GetEventFromBuffer(char* ap_buffer, int a_th)
{
  #ifdef CASTOR_VERBOSE
  // Verbose
  if(m_verbose>=4) Cout("iDataFileSPECT::GetEventFromBuffer() ..." << endl);
  #endif

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

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckSpecificParameters()
  \brief  Check parameters specific to SPECT data
  \return 0 if success, and positive value otherwise.
*/
int iDataFileSPECT::CheckSpecificParameters()
{
  // Verbose
  if (m_verbose>=3) Cout("iDataFileSPECT::CheckSpecificParameters()" << endl);
  
  // Error if m_dataType != SPECT
  if (m_dataType != TYPE_SPECT)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Error, Data type should be SPECT (=1) !'" << endl);
    return 1;
  }

  if (m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Error, number of projections not initialized (should be >0) !'" << endl);
    return 1;
  }
  
  if (mp_angles == NULL)
  {
    Cerr("***** iDataFileSPECT::CheckSpecificParameters() -> Error, Projection angles not initialized !'" << endl);
    return 1;
  }
  
  
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckFileSizeConsistency()
  \brief  This function is implemented in child classes
          Check if file size is consistent.
  \return 0 if success, and positive value otherwise.
*/
int iDataFileSPECT::CheckFileSizeConsistency()
{
  if(m_verbose >=3) Cout("iDataFileSPECT::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("***** iDataFileSPECT::CheckFileSizeConsistency() -> Error : Datafile size is not consistent with the information provided by the user/datafile!" << endl);
    Cerr("***** iDataFileSPECT::CheckFileSizeConsistency() -> Expected size : "<< m_totalNbEvents*m_sizeEvent << endl);
    Cerr("***** iDataFileSPECT::CheckFileSizeConsistency() -> Actual size : "<< sizeInBytes << endl << endl);

    Cerr("***** iDataFileSPECT::CheckFileSizeConsistency() -> ADDITIONAL INFORMATION ABOUT THE DATAFILE INITIALIZATION : " << endl);


    if(m_eventKindFlag)
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> Coincidence kind (trues, scatter, ...) term is enabled " << endl);
    else
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> No informations about the kind of coincidences in the data " << endl);
      
    if(m_normCorrectionFlag)
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> Normalization correction term is enabled " << endl);
    else
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> No normalization correction terms in the data " << endl);
      
    if(m_scatCorrectionFlag)
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> Scatter correction term is enabled " << endl);
    else
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> No scatter correction terms in the data " << endl);

    Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> Calibration factor value is : " << m_calibrationFactor << endl);
    Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> Number of bins are equal to "<<mp_nbOfBins[0]<<","<< mp_nbOfBins[1]<<" for X,Y axis respectively." << endl);

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

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

    if(mp_POIResolution[2]<0)
      Cerr("      iDataFileSPECT::CheckFileSizeConsistency() -> No POI enabled on the Z axis " << endl);
    else
      Cerr("      iDataFileSPECT::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_eventKindFlag)
      Cout("  --> Coincidence kind (trues, scatter, ...) 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_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 X,Y axis respectively." << 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;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn      iDataFileSPECT::InitAngles()
  \param   ap_angles
  \brief   allocate memory for the mp_angles variable using m_nbProjections
           and initialize the projection angles with the provided list of values
  \return  0 if success, positive value otherwise
*/
int iDataFileSPECT::InitAngles(FLTNB* ap_angles)
{
  if(m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::InitAngles() -> Error, number of projection angles not initialized !'" << endl);
    return 1;
  }
  
  mp_angles = new FLTNB[m_nbOfProjections];
  
  for(size_t a=0 ; a<m_nbOfProjections ; a++)
    mp_angles[a] = ap_angles[a];

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn      iDataFileSPECT::InitCorToDetectorDistance()
  \param   ap_CORtoDetectorDistance
  \brief   allocate memory for the ap_CORtoDetectorDistance variable
           using m_nbProjections, and initialize the projection angles
           with the provided list of values
  \return  0 if success, positive value otherwise
*/
int iDataFileSPECT::InitCorToDetectorDistance(FLTNB* ap_CORtoDetectorDistance)
{
  if(m_nbOfProjections == 0)
  {
    Cerr("***** iDataFileSPECT::InitAngles() -> Error, number of projection angles not initialized !'" << endl);
    return 1;
  }

  mp_CORtoDetectorDistance = new FLTNB[m_nbOfProjections];

  for(size_t a=0 ; a<m_nbOfProjections ; a++)
    mp_CORtoDetectorDistance[a] = ap_CORtoDetectorDistance[a];

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