source: Sophya/trunk/AddOn/TAcq/brproc.cc@ 4065

//----------------------------------------------------------------
// Projet BAORadio - (C) LAL/IRFU  2008-2011
// Classes de threads de traitement pour BAORadio 
//----------------------------------------------------------------

#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <fstream>
#include <signal.h>

#include "pexceptions.h"
#include "tvector.h"
#include "ntuple.h"
#include "datatable.h"
#include "histos.h"
#include "fioarr.h"
#include "matharr.h"
#include "timestamp.h"
#include "ctimer.h"
#include "fftpserver.h"
#include "fitsarrhand.h"

#include "FFTW/fftw3.h"


#include "pciewrap.h"
#include "brpaqu.h"
#include "brproc.h"



//---------------------------------------------------------------------
// Classe BRMeanSpecCalculator de traitement de spectres - 
// Calcul de spectres moyennes,variance / voie + nettoyage
// Implementation de traitement par fenetres temps-frequence
// Le temps correspond au numero de paquet
//---------------------------------------------------------------------
/*!
   \class BRMeanSpecCalculator
   \ingroup TAcq

   \brief Computes mean spectra on FFT data.
*/

/* --Methode-- */
BRMeanSpecCalculator::BRMeanSpecCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, 
                                           bool fgdatafft, bool fgsinglechan)
  : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean),
    fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan), 
    nbpaq4mean_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()), 
    nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres())
{
  setNameId("meanSpecCalc",1);

  uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
  const char* venvp = NULL; 
  venvp=getenv("BRANA_NBYTECUT");
  if (venvp!=NULL){
    nb_octets_entrop = atoi(venvp);
    cout << "BRMeanSpecCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
  }

  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);

  if (fgsinglechannel_) {
    mspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); 
    sigspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); 
    sgain_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); 
  }
  else {
    mspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); 
    sigspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); 
    sgain_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); 
  }
  mspecmtx_=(r_4)(0.);
  sigspecmtx_=(r_4)(0.);
  sgain_=(r_4)(1.);    // Gain en fonction de la frequence, à 1 par defaut 

  numfile_=0;
  totnbpaq_=0;
  moyfc_=moytt_=0.;
  nbpmoyttfc_=0;

  size_t nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
  
  for(size_t i=0; i<nchan; i++) {
    nbpaq4mean_[i]=nbadpaq_[i]=0;
  }
  

  // Definition des tailles de fenetres de spectres, etc ...
  SetSpectraWindowSize();
  SetMaxNbSpecWinFiles();
  nbtot_specwin_=0;
  SetVarianceLimits();
  SetNumberOfBands();

  ofsdtp_=NULL;
  dtp_=NULL;
  xnt_=NULL;
  ofsdtpms_=NULL;
  dtpms_=NULL;
  xntms_=NULL;
  firstfreqms_=0; 
  lastfreqms_=-1;
  nbandms_=-1;
}

/* --Methode-- */
BRMeanSpecCalculator::~BRMeanSpecCalculator()
{
  cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
  uint_8 npqm=0;
  for(size_t i=0; i<nbpaq4mean_.size(); i++)  npqm+=nbpaq4mean_[i];
  if (npqm>nmean_*nbpaq4mean_.size()/10)  SaveMeanSpectra();
  for(size_t i=0; i<nbadpaq_.size(); i++) {
    cout << " Channel " << i << "  NBadPaq=" << nbadpaq_[i] << " / TotNbPaq=" << totnbpaq_ << endl;
  }
  if (dtp_) {
    cout << *dtp_; 
    delete dtp_;
    delete ofsdtp_;
  }
  if (dtpms_) {
    cout << *dtpms_; 
    delete dtpms_;
    delete ofsdtpms_;
  }
  if (xnt_)   delete xnt_;
  if (xntms_)   delete xntms_;
  cout << " ------------------------------------------------------------------------ " << endl;
}

/* --Methode-- */
void BRMeanSpecCalculator::SetSpectraWindowSize(uint_4 winsz, uint_4 wszext)
{
  if (winsz < 3) {
    winsz=1;  wszext=0;
  }  
  if (wszext>=winsz/2) wszext=winsz/2;
  sa_size_t sz[5]={0,0,0,0,0};
  sz[0]=mspecmtx_.NCols();
  sz[1]=mspecmtx_.NRows();
  sz[2]=winsz+2*wszext;
  spec_window_.SetSize(3,sz);
  spwin_ext_sz_=wszext;
  sz[0]=mspecmtx_.NRows();
  sz[1]=winsz+2*wszext;  
  clnflg_.SetSize(2,sz);
  cout << "BRMeanSpecCalculator::SetSpectraWindowSize()/Info: SpectraWindowSize()=" << GetSpectraWindowSize()
       << " ExtensionSize=" << GetSpecWinExtensionSize() << " Overlap=" << GetSpecWinOverlapSize() 
       << " ArraySize=" << spec_window_.SizeZ() << endl;

  paqnum_w_start=spwin_ext_sz_;  // premiere initialisation du numero de paquet 
  return;
}

/* --Methode-- */
void BRMeanSpecCalculator::DefinePerPaquetDataTable()
{
  cout << "BRMeanSpecCalculator::DefinePerPaquetDataTable: Creating DataTable" << endl;
  string dtfile="!"+outpath_+"/dtspecpaq.fits";
  ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
  dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true);
  char cnom[32];
  size_t nchan=mspecmtx_.NRows();
  for(int i=0; i<nchan; i++) {
    sprintf(cnom,"variance%d",i);
    dtp_->AddFloatColumn(cnom);    
  }
  for(int i=0; i<nchan; i++) {
    for(int j=0;j<numberOfBands_;j++){
      sprintf(cnom,"varnormb%d%d",i,j);
      dtp_->AddFloatColumn(cnom);    
    }
  }
  /*
    for(int i=0; i<nchan; i++) {
    sprintf(cnom,"sigma%d",i);
    dtp_->AddFloatColumn(cnom);    
    }
  */
  //  xnt_=new double[nchan*2];  CHECK : faut-il reallouer ?
  cout << " ...Number of columns " << nchan*(1+numberOfBands_) << " doubles -> to FITS file" << dtfile << endl;
  xnt_=new double[nchan*(1+numberOfBands_)];

}

/* --Methode-- */
void BRMeanSpecCalculator::DefineTimeAvgPowerDataTable(int first, int last, int nband)
{
  cout << "BRMeanSpecCalculator::DefineTimeAvgPowerDataTable: Creating time averged power spectrum" << endl;
  string dtfile="!"+outpath_+"/avgspec.fits";
  ofsdtpms_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
  dtpms_ = new SwFitsDataTable(*ofsdtpms_,1024,true);
  dtpms_->AddDoubleColumn("mtimetag");
  dtpms_->AddDoubleColumn("time");
  
  if (last>first) {
    if ((first<0)||(first>=mspecmtx_.NCols()))  first=0;
    if ((last<0)||(last>mspecmtx_.NCols()))  last=mspecmtx_.NCols();
    if (nband<1)  nband=1;
  }
  else {
    first=1;  last=mspecmtx_.NCols();  nband=1;
  }
  firstfreqms_=first;   lastfreqms_=last;  nbandms_=nband;
  char cnom[32];
  size_t nchan=mspecmtx_.NRows();
  for(int i=0; i<nchan; i++) {
    for(int j=0; j<nband; j++) {
      sprintf(cnom,"pow%dc%d",j,i);
      dtpms_->AddFloatColumn(cnom); 
    }   
  }
  cout << " ...Number of columns " << (2+nchan*nband) << " doubles -> to FITS file" << dtfile << endl;
  xntms_=new double[nchan*nband+2];
}
  
/* --Methode-- */
void BRMeanSpecCalculator::SetNumberOfBands(int numberOfBands, int ibandfirst, int ibandlast)
{
  if (numberOfBands < 1) numberOfBands = 1;
  numberOfBands_ = numberOfBands;
  if (ibandfirst < 0 )ibandfirst = 0;
  if (ibandlast >= numberOfBands_ ) ibandlast = numberOfBands_-1;
  ibandfirst_=ibandfirst; ibandlast_=ibandlast;

  cout << "BRMeanSpecCalculator::SetNumberOfBands (END) : "
       << numberOfBands_ << " "
       << ibandfirst_ << " "
       << ibandlast_ << endl;

}

/* --Methode-- */
void BRMeanSpecCalculator::ReadGainFitsFile(string filename, bool fgapp)
{
  cout << " BRMeanSpecCalculator::ReadGainFitsFile() - reading file " << filename;
  FitsInOutFile fis(filename, FitsInOutFile::Fits_RO);
  fis >> sgain_;
  fg_apply_gains_=fgapp;
  cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << " ApplyGains=" 
       << ((fg_apply_gains_)?"true":"false") << endl;
  if( (spec_window_.SizeX()!= sgain_.NCols()) || (spec_window_.SizeY()!= sgain_.NRows()) ){
    cout << " ReadGainFitsFile: BAD Gain Matrix sizes " << endl;
    sgain_.Show();
    spec_window_.Show();
    throw ParmError("ReadGainFitsFile: BAD Gain Matrix sizes");
  }
}

//JEC
//static inline r_4 Zmod2(complex<r_4> z) 
//{ return (z.real()*z.real()+z.imag()*z.imag()); }




/* --Methode-- */
int BRMeanSpecCalculator::Process()
{
  // Cette methode remplit le tableau spec_window_ avec les spectres (renormalise avec 
  // les gains si demande) et appelle la methode du traitement de la fenetre temporelle 
  // des spectres le cas echeant ProcSpecWin()

 
  int_8 nbpaqdec = (int_8)totnbpaq_-(int_8)GetSpecWinOverlapSize();
  if ((nbpaqdec>0)&&(nbpaqdec%GetSpectraWindowSize()==0)) {
    paqnum_w_end=totnbpaq_-GetSpecWinExtensionSize();
    ProcSpecWin(paqnum_w_start, paqnum_w_end);
    paqnum_w_start=totnbpaq_-GetSpecWinExtensionSize();
  }

  if (fgdatafft_) {  // Donnees firmware FFT 
    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
      TwoByteComplex* zp=NULL;
      if (fgsinglechannel_) {
        zp=vpaq_[i].Data1C();
      }
      else {
        zp=vpaq_[i/2].Data1C();
        if (i%2==1)  zp=vpaq_[i/2].Data2C();
      }
      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
      for(sa_size_t f=0; f<spec_window_.SizeX(); f++) 
        spec_window_(f,i,kz) = zp[f].module2F();
    }
  }
  else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator
    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
      complex<ODT>* zp=NULL;
      if (fgsinglechannel_) {
        zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i]);
      }
      else {
        zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i/2]);
        if (i%2==1)  zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ;
      }
      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
      for(sa_size_t f=0; f<spec_window_.SizeX(); f++) 
        spec_window_(f,i,kz) = Zmod2(zp[f]); 
    }
  }
  if (fg_apply_gains_) {   // Application des gains, si demande 
    sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
      (spec_window_(Range::all(), Range(i), Range(kz)).CompactAllDimensions()).Div(sgain_.Row(i).CompactAllDimensions());
    }
  }

  // Pour calcul MeanTimeTag , MeanFrameCounter 
  moyfc_+=curfc_[0];
  moytt_+=((double)vpaq_[0].TimeTag()/1.25e8);
  nbpmoyttfc_++;
  totnbpaq_++;
  return 0;
}


/* --Methode-- */
void BRMeanSpecCalculator::ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend)
{

  if (prtlev_>0)  {
    uint_8 modulo =  prtmodulo_/GetSpectraWindowSize();
    if (modulo<1) modulo=1;
    if (nbtot_specwin_%modulo==0) {
      cout << " BRMeanSpecCalculator::ProcSpecWin() num_win=" << nbtot_specwin_ << " numpaqstart=" << numpaqstart 
           << " numpaqend=" << numpaqend << endl;
      cout << " ... ObsTime=" << getObsTime() << " TimeTag=" << getCurTimeTagSeconds() << " s. FrameCounter=" 
           << getCurFrameCounter() << endl;
    }
  }
  // On appelle la routine de nettoyage qui doit flagger les mauvais paquets 
  FlagBadPackets(numpaqstart, numpaqend);

  // Boucle sur les numeros de paquets de la fenetre en temps
  for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
    // On sauvegarde les spectres moyennes si necessaire 
    if ((nbpaq4mean_[0]>0)&&(nbpaq4mean_[0]%nmean_ == 0))  SaveMeanSpectra();
    //  On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
    sa_size_t kz=PaqNumToArrayIndex(jp); 
    // Boucle sur les numeros de voie (canaux) 
    for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  {
      if ( clnflg_(i,kz) != 0)  continue;
      TVector< r_4 > spec = mspecmtx_.Row(i); 
      TVector< r_4 > sspec = sigspecmtx_.Row(i); 
      // Calcul de spectres moyennes et variance 
      for(sa_size_t f=1; f<spec.Size(); f++)  {   // boucle sur les frequences 
        spec(f)  += spec_window_(f,i,kz);
        sspec(f) += spec_window_(f,i,kz)*spec_window_(f,i,kz);  
      }
      nbpaq4mean_[i]++;  // compteur de paquets OK pour la moyenne 
    }
  }
  if (nbtot_specwin_<nmaxfiles_specw_)  SaveSpectraWindow();
  nbtot_specwin_++;
  return;
}

/* --Methode-- */
void BRMeanSpecCalculator::FlagBadPackets(uint_8 numpaqstart, uint_8 numpaqend)
{
  // Boucle sur les numeros de paquets de la fenetre en temps
  for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
    //  On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
    sa_size_t kz=PaqNumToArrayIndex(jp); 
    // Boucle sur les numeros de voie (canaux) 
    for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  {

      clnflg_(i,kz)=0;

      double mean, sigma;
      ////////BUG      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
      double variance=0.;
      variance=spec_window_(Range(1,Range::lastIndex()), Range(i), Range(kz)).Sum();
      if(xnt_)xnt_[i]=variance;
      if(numberOfBands_>0){
        //Compute nomalized variance in bands freq. 
        sa_size_t fMin;
        sa_size_t fMax;
        int bandW = spec_window_.SizeX()/numberOfBands_;
        vector<double> varNomBinned(numberOfBands_);
        for (sa_size_t j=ibandfirst_; j<=ibandlast_; j++){
          fMin = j*bandW;
          fMax =fMin+bandW-1;
          varNomBinned[j]=spec_window_(Range(fMin,fMax), Range(i), Range(kz)).Sum();
          varNomBinned[j]/=(r_4)bandW;
          if(xnt_)xnt_[spec_window_.SizeY()+i*numberOfBands_+j] =  varNomBinned[j];

          //        cout << "(jec) var["<<j<<"] =" << varNomBinned[j]
          //     << " min " << varmin_
          //     << " max " << varmax_ << endl;
          if(varNomBinned[j]<varmin_)      { clnflg_(i,kz)=10+j;  nbadpaq_[i]++; break;}
          else if(varNomBinned[j]>varmax_) { clnflg_(i,kz)=100+j; nbadpaq_[i]++; break;}
        }
        //cout << "clnflg_("<<i<<","<<kz<<"): " << clnflg_(i,kz)  << endl;
      }//if bands
    }//loop on channels

    if (dtp_ && xnt_)  dtp_->AddRow(xnt_);
  }
  return;
}

/* --Methode-- */
void BRMeanSpecCalculator::SaveMeanSpectra()
{
  for(sa_size_t ir=0; ir<mspecmtx_.NRows(); ir++) {
    char buff[32];
    sprintf(buff,"NPAQSUM_%d",(int)ir);
    mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
    mspecmtx_.Info()[buff] = nbpaq4mean_[ir];
    sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
    sigspecmtx_.Info()[buff] = nbpaq4mean_[ir];
    if (nbpaq4mean_[ir] > 0) {
      mspecmtx_.Row(ir)  /= (r_4)nbpaq4mean_[ir];
      sigspecmtx_.Row(ir) /= (r_4)nbpaq4mean_[ir];
      sigspecmtx_.Row(ir) -= (mspecmtx_.Row(ir) && mspecmtx_.Row(ir));  // Mean(X^2) - [ Mean(X) ]^2
    }
  }
  mspecmtx_.Info()["ENDTIME"]=getObsTime();
  sigspecmtx_.Info()["ENDTIME"]=getObsTime();
  if (nbpmoyttfc_>1)  {
    moyfc_/=(double)nbpmoyttfc_;
    moytt_/=(double)nbpmoyttfc_;
  }
  string ikey,ikdesc;    
  ikey="MeanFC";  ikdesc="Mean FrameCounter";
  mspecmtx_.Info().SetI(ikey,moyfc_); 
  mspecmtx_.Info().SetComment(ikey,ikdesc);
  sigspecmtx_.Info().SetI(ikey,moyfc_); 
  sigspecmtx_.Info().SetComment(ikey,ikdesc);
  
  ikey="MeanTT";  ikdesc="Mean TimeTag";
  mspecmtx_.Info().SetD(ikey,moytt_);
  mspecmtx_.Info().SetComment(ikey,ikdesc);
  sigspecmtx_.Info().SetD(ikey,moytt_);
  sigspecmtx_.Info().SetComment(ikey,ikdesc);
  
  char nfile[64];
  string flnm;
  {
  sprintf(nfile,"mspecmtx%d.fits",numfile_);
  flnm="!"+outpath_+nfile;
  FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
  fos << mspecmtx_;
  }
  {
  sprintf(nfile,"sigspecmtx%d.fits",numfile_);
  flnm="!"+outpath_+nfile;
  FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
  fos << sigspecmtx_;
  }

  cout << numfile_ << "-BRMeanSpecCalculator::SaveMeanSpectra() NPaqProc="
       << totnbpaq_ << "  -> Mean/Sig spectra Matrix in " << flnm << " /sigspec...ppf" << endl;
  numfile_++;  

  if (dtpms_!=NULL)  FillPwrTmDTable(mspecmtx_);

  for(size_t i=0; i<nbpaq4mean_.size(); i++) nbpaq4mean_[i]=0;
  mspecmtx_ = (r_4)(0.);
  sigspecmtx_ = (r_4)(0.);
  moyfc_=moytt_=0.;  
  nbpmoyttfc_=0;
  return;
}

/* --Methode-- */
void BRMeanSpecCalculator::SaveSpectraWindow()
{
  char nfile[64];
  string flnm;
  sprintf(nfile,"specwin%d.fits",nbtot_specwin_);
  flnm="!"+outpath_+nfile;
  FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
  fos << spec_window_;
  cout << " SaveSpectraWindow() " << nbtot_specwin_ << "- file " << nfile << " created " << endl;
}

/* --Methode-- */
void BRMeanSpecCalculator::FillPwrTmDTable(TMatrix< r_4 >& specmtx)
{
  if (dtpms_==NULL) return;
  xntms_[0]=moytt_;  
  xntms_[1]=getObsTimeSeconds();
  size_t knt=2;  
  sa_size_t wband=(lastfreqms_-firstfreqms_)/nbandms_;
  for(sa_size_t ir=0; ir<specmtx.NRows(); ir++) {
    for(sa_size_t jf=0; jf<specmtx.NCols()-wband; jf+=wband) {
      TVector<r_4> srow=specmtx.Row(ir);
      xntms_[knt]=srow.SubVector(Range(jf,jf+wband-1)).Sum();
      knt++;
    }
  }
  dtpms_->AddRow(xntms_);
}

//---------------------------------------------------------------------
// Classe de thread de calcul de FFT sur donnees RAW 
//---------------------------------------------------------------------
/*!
   \class BRFFTCalculator
   \ingroup TAcq

   \brief Computes FFT on data from RAW ADC board firmware.
*/
/* --Methode-- */
BRFFTCalculator::BRFFTCalculator(RAcqMemZoneMgr& memgr, bool fgsinglechannel)
  : BRBaseProcessor(memgr), fgsinglechannel_(fgsinglechannel), totnbfftpaq_(0)
{
  uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
  const char* venvp = NULL; 
  venvp=getenv("BRANA_NBYTECUT");
  if (venvp!=NULL){
    nb_octets_entrop = atoi(venvp);
    cout << "BRFFTCalculator  : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
  }

  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
  //JEC END
  //  BRPaquet paq(memgr_.PaqSize());
  setNameId("FFTCalc",2);
  ffts_.SetInDataSize((fgsinglechannel_)?paq.DataSize():paq.DataSize()/2);
}

/* --Methode-- */
BRFFTCalculator::~BRFFTCalculator()
{
}


/* --Methode-- */
int BRFFTCalculator::Process()
{
  for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
    ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data1() ),  
                 reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]) );
    totnbfftpaq_++;
    if ( fgsinglechannel_ ) continue;
    ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data2() ),  
                 reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]+memgr_.ProcPaqSize()/2) );
    totnbfftpaq_++;
  }    
  return 0;
}


//-------------------------------------------------------------------------
// Classe WBRFFT : Calcul de TF sur donnees brutes (firmware RAW)
//-------------------------------------------------------------------------
/*!
   \class WBRFFT
   \ingroup TAcq

   \brief Helper class for FFT computation used by BRFFTCalculator class.

   This class uses FFTW library.
*/

ZMutex* WBRFFT::mtx_fftwp_ = NULL;

/* --Methode-- */
WBRFFT::WBRFFT(uint_4 sz)
  : sz_(sz)
{
  if (mtx_fftwp_ == NULL) mtx_fftwp_ = new ZMutex;
  if (sz>0)  SetInDataSize(sz);
}

/* --Methode-- */
WBRFFT::~WBRFFT()
{
}

/* --Methode-- */
void WBRFFT::SetInDataSize(uint_4 sz)
{
  sz_ = sz; 
  if (sz_<1) return;
  inp.SetSize(sz);
  outfc.SetSize(sz/2+1);
  mtx_fftwp_->lock();
  myplan_ = fftwf_plan_dft_r2c_1d(inp.Size(), inp.Data(), 
                       (fftwf_complex*)outfc.Data(), FFTW_ESTIMATE); 
  mtx_fftwp_->unlock();
}

/* --Methode-- */
void WBRFFT::DoFFT( IDT *indata, complex<ODT> * ofc)
{
  if (sz_<1) return;
  for(uint_4 k=0; k<inp.Size(); k++)   inp(k)=(ODT)indata[k];
  fftwf_execute(myplan_); 
  for(uint_4 k=0; k<outfc.Size(); k++)   ofc[k]=outfc(k)/(ODT)sz_;   // on renormalise les coeff FFT ( / sz )
  return;
}

/* --Methode-- */
void WBRFFT::PrintData(IDT *indata,  complex<ODT> * ofc, uint_4 sz)
{
    cout << " --- WBRFFT::PrintData() size=" << sz << endl;
  for(uint_4 k=0; k<sz; k+=8) {
    IDT* in = indata+k;
    cout << " Indata[" << k << "..." << k+8 << "]= ";
    for(uint_4 i=0; i<8; i++)   cout << (IIDT)in[i] << "  ";
    cout << endl;  
  }
  cout << endl;
  for(uint_4 k=0; k<sz/2; k+=4) {
    complex< ODT>* out = ofc+k;
    cout << " OutFC[" << k << "..." << k+4 << "]= ";
    for(uint_4 i=0; i<4; i++)   cout << out[i] << "  ";
    cout << endl;  
  }
  return;

}


//---------------------------------------------------------------
// Classe thread de traitement donnees ADC avec 2 voies par frame
//         !!!! OBSOLETE  !!!!   
//---------------------------------------------------------------
/*!
   \class BRProcA2C
   \ingroup TAcq

   \brief Deprecated processing class. 
*/

// Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe
static ZMutex* pmutfftw=NULL;

/* --Methode-- */
BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean, 
                           uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card)
  :  memgr(mem) 
{
  fgraw_ = fgraw;
  nmax_ = nmax; 
  nmean_ = nmean;
  if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl;
  else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl;
  nfsmap_ = nfsmap;
  stop_ = false;        
  path_ = path; 
  fgnotrl_ = fgnotrl;
  fghist_ = fghist;
  card_ = card;
  if (pmutfftw==NULL) pmutfftw=new ZMutex;  
}

/* --Methode-- */
void BRProcA2C::Stop()
{
 stop_=true;
 // cout <<" BRProcA2C::Stop ... > STOP " << endl;
}




static inline string card2name_(int card)
{
  if (card==2) return " (Chan3,4) ";
  else return " (Chan1,2) ";
}
/* --Methode-- */
void BRProcA2C::run()
{
  setRC(1);     
  try {
    Timer tm("BRProcA2C", false);
    TimeStamp ts; 
    BRPaqChecker pcheck(!fgnotrl_);  // Verification/comptage des paquets 
    
    size_t totnbytesout = 0;
    size_t totnbytesproc = 0;

    cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ 
         << " NMean=" << nmean_ << card2name_(card_) << endl;   
    cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl;
    char fname[512];
//    sprintf(fname,"%s/proc.log",path_.c_str());
//    ofstream filog(fname);
//    filog << " BRProcA2C::run() - starting log file " << ts << endl;                 
//    filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;       

/*----DELETE   NTuple 
    const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"};
    NTuple nt(8, nnames);
    double xnt[10];
    uint_4 nmnt = 0;
    double ms1,ms2,ms12,ms12re,ms12im,ms12phi;
----*/
// Time sample (raw data) /FFT coeff histograms
   Histo* ph1=NULL;
   Histo* ph2=NULL;
   if (fghist_) {
     if (fgraw_) {
       ph1 = new Histo(-0.5, 255.5, 256);     
       ph2 = new Histo(-0.5, 255.5, 256);     
     }
     else {
       ph1 = new Histo(-128.5, 128.5, 257);     
       ph2 = new Histo(-128.5, 128.5, 257);     
     }
   }

// Initialisation pour calcul FFT 
    TVector< complex<r_4> > cfour1;  // composant TF
    uint_4 paqsz = memgr.PaqSize();
    uint_4 procpaqsz = memgr.ProcPaqSize();
    
   
    BRPaquet pq(NULL, NULL, paqsz); 
    TVector<r_4> vx(pq.DataSize()/2);
    int szfour = pq.DataSize()/2/2+1;
    cfour1.SetSize(szfour);
/*
    vx = (r_4)(0.);
    FFTPackServer ffts;
    ffts.FFTForward(vx, cfour1);
    szfour = cfour1.Size();
*/

    bool fgtimfreq = false;  // true->cartes temps<>frequences
    if (nfsmap_>0) fgtimfreq=true;

    TVector< complex<r_4> > cfour2(cfour1.Size());
    
    TVector<r_4> spectreV1(cfour1.Size());
    TVector<r_4> spectreV2(cfour1.Size());
    TVector<r_4> moyspecV1(cfour1.Size());   // Moyenne des Spectres 
    TVector<r_4> moyspecV2(cfour1.Size());
    TVector<r_4> sigspecV1(cfour1.Size());   // Sigma des Spectres
    TVector<r_4> sigspecV2(cfour1.Size());
    TVector< complex<r_4> > visiV12( cfour1.Size() );

    TMatrix<r_4> timfreqV1, timfreqV2;   // Cartes temps<>frequences
    if (fgtimfreq) {
      timfreqV1.SetSize(nmean_, spectreV1.Size()/nfsmap_);
      timfreqV2.SetSize(nmean_, spectreV2.Size()/nfsmap_);
    }
    cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz 
         << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
         << " *8="  << cfour1.Size()*8 << endl;

    pmutfftw->lock();
    fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), 
                          (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE); 
    fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), 
                          (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE); 
    pmutfftw->unlock();

    uint_4 ifile = 0;                           
    uint_4 nzm = 0;  // Nb de paquets moyennes pour le calcul de chaque spectre
    uint_4 nmoyspec = 0;  // Nb de spectres moyennes

    uint_4 curfc=0;
    uint_8 curtt=0;
    uint_8 firsttt=0;
    bool fgfirst=true;
    double moysig[2]={0.,0.};
    double sigsig[2]={0.,0.};
    uint_8 nbsig[2]={0,0};

    for (uint_4 kmz=0; kmz<nmax_; kmz++) {
      if (stop_) break;
      int mid = memgr.FindMemZoneId(MemZA_ProcA);
      Byte* buff = memgr.GetMemZone(mid);
      if (buff == NULL) {
         cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
             break;             
      }
      Byte* procbuff = memgr.GetProcMemZone(mid);
      if (procbuff == NULL) {
            cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
        break;  
      }
//---- DELETE      nmnt=0;  ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
      for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
            BRPaquet paq(NULL, buff+i*paqsz, paqsz); 
        if (!pcheck.Check(paq)) continue;   // on ne traite que les paquets OK
        if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; } 
        curfc=paq.FrameCounter();
        curtt=paq.TimeTag()-firsttt;
// Traitement voie 1        
        if (fghist_) {
          for(sa_size_t j=0; j<vx.Size(); j++) {
            r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j)));
            ph1->Add((r_8)vts);
            moysig[0] += (double)vts;
            sigsig[0] += ((double)vts)*((double)vts);
            nbsig[0]++;
          }
          for(sa_size_t j=0; j<vx.Size(); j++) {
            r_4 vts=(fgraw_)?((r_4)(*(paq.Data2()+j))):((r_4)(*(paq.Data2S()+j)));
            ph2->Add((r_8)vts);
            moysig[1] += (double)vts;
            sigsig[1] += ((double)vts)*((double)vts);
            nbsig[1]++;
          }
        }
        if (fgraw_) {
          for(sa_size_t j=0; j<vx.Size(); j++) 
            vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
          //        fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz);
          fftwf_execute(plan1); 
          // Traitement voie 2              
          for(sa_size_t j=0; j<vx.Size(); j++) 
            vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
          fftwf_execute(plan2); 
        }
        else {
          for(sa_size_t j=1; j<cfour1.Size()-1; j++) {
            cfour1(j) = complex<r_4>((r_4)paq.Data1C()[j].realB(), (r_4)paq.Data1C()[j].imagB());
            cfour2(j) = complex<r_4>((r_4)paq.Data2C()[j].realB(), (r_4)paq.Data2C()[j].imagB());
          }
          cfour1(0) = complex<r_4>((r_4)paq.Data1C()[0].realB(), (r_4)0.);
          cfour1(cfour1.Size()-1) = complex<r_4>((r_4)paq.Data1C()[0].imagB(), (r_4)0.);
          cfour2(0) = complex<r_4>((r_4)paq.Data2C()[0].realB(), (r_4)0.);
          cfour2(cfour2.Size()-1) = complex<r_4>((r_4)paq.Data2C()[0].imagB(), (r_4)0.);
        }
        for(sa_size_t j=0; j<spectreV1.Size(); j++) 
          spectreV1(j) += BRMeanSpecCalculator::Zmod2(cfour1(j));
        memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
        if (fgtimfreq) {   // Remplissage tableau temps-frequence
          for(sa_size_t c=1; c<timfreqV1.NCols(); c++) {
            for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) 
              timfreqV1(nzm, c) += BRMeanSpecCalculator::Zmod2(cfour1(j));
          }
        }
        for(sa_size_t j=0; j<spectreV2.Size(); j++) 
          spectreV2(j) += BRMeanSpecCalculator::Zmod2(cfour2(j));  // BRMeanSpecCalculator::Zmod2(zp2[j]); 
        memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
        if (fgtimfreq) {   // Remplissage tableau temps-frequence
          for(sa_size_t c=1; c<timfreqV2.NCols(); c++) {
            for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) 
              timfreqV2(nzm,c) += BRMeanSpecCalculator::Zmod2(cfour2(j));
          }
        }

// Calcul correlation (visibilite V1 * V2)
        for(sa_size_t j=0; j<visiV12.Size(); j++) 
          visiV12(j)+=cfour1(j)*conj(cfour2(j));
//        for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j];
        if (nzm==0) {
          spectreV1.Info()["StartFC"] = curfc;
          spectreV2.Info()["StartFC"] = curfc;
          visiV12.Info()["StartFC"] = curfc;
          spectreV1.Info()["StartTT"] = curtt;
          spectreV2.Info()["StartTT"] = curtt;
          visiV12.Info()["StartTT"] = curtt;
        }
        nzm++;   
/*----DELETE
        if (nmnt==0)  { xnt[0]=paq.FrameCounter();  xnt[1]=paq.TimeTag(); }
        for(sa_size_t j=2700; j<2800; j++) { 
          ms1 += BRMeanSpecCalculator::Zmod2(cfour1(j)); ms2 += BRMeanSpecCalculator::Zmod2(cfour2(j)); 
          complex<r_4> zvis =  cfour1(j)*conj(cfour2(j));
          ms12 += BRMeanSpecCalculator::Zmod2(zvis);   ms12re += zvis.real();  ms12im += zvis.imag();
          ms12phi+= atan2(zvis.imag(),zvis.real());
        } 
        nmnt++;    
----*/
        totnbytesproc += paq.DataSize();
        totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());

      } // Fin de boucle sur les paquets d'une zone

/*---- DELETE
      if (nmnt>0)  {
        double fnorm = (double)nmnt*(2800-2700); 
        xnt[2] = ms1 /= fnorm;
        xnt[3] = ms2 /= fnorm;
        xnt[4] = ms12 /= fnorm;
        xnt[5] = ms12re /= fnorm;
        xnt[6] = ms12im /= fnorm;
        xnt[7] = ms12phi /= fnorm;
        nt.Fill(xnt);
      }
----*/
      if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
        spectreV1 /= (r_4)(nzm);
        spectreV2 /= (r_4)(nzm);

        // pour le calcul des moyennes et sigmas de ces spectres 
        moyspecV1 += spectreV1;
        moyspecV2 += spectreV2;
        sigspecV1 += (spectreV1 && spectreV1);
        sigspecV2 += (spectreV2 && spectreV2);
        nmoyspec++;

        visiV12 /= complex<r_4>((r_4)nzm, 0.);

        spectreV1.Info()["NPaqMoy"] = nzm;
        spectreV2.Info()["NPaqMoy"] = nzm;
        visiV12.Info()["NPaqMoy"] = nzm;
        spectreV1.Info()["EndFC"] = curfc;
        spectreV2.Info()["EndFC"] = curfc;
        visiV12.Info()["EndFC"] = curfc;
        spectreV1.Info()["EndTT"] = curtt;
        spectreV2.Info()["EndTT"] = curtt;
        visiV12.Info()["EndTT"] = curtt;
        {
        sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
        POutPersist po(fname);
        string tag1="specV1";
        string tag2="specV2";
        string tag12="visiV12";
        string tagh1="tshV1";
        string tagh2="tshV2";
        string tagtf1="timfreqV1";
        string tagtf2="timfreqV2";
        if (card_==2) {
          tag1 = "specV3";
          tag2 = "specV4";
          tagh1 = "tshV1";
          tagh2 = "tshV2";
          tag12="visiV34";
          tagtf1="timfreqV3";
          tagtf2="timfreqV4";
        }
        po << PPFNameTag(tag1) << spectreV1; 
        po << PPFNameTag(tag2) << spectreV2; 
        po << PPFNameTag(tag12) << visiV12; 
        if (fghist_) {
          po << PPFNameTag(tagh1) << (*ph1); 
          po << PPFNameTag(tagh2) << (*ph2); 

          double sspvmax[3] = {0.,0.,0.};
          int_4 sspvmaxidx[3] = {-1,-1,-1};
          for(int jji=1;jji<visiV12.Size()-1;jji++) {
            r_4 zmv2 = BRMeanSpecCalculator::Zmod2(visiV12(jji)); 
            if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; }
          }
          TVector<r_4>& sspv = spectreV1;
          for(int ic=0; ic<2; ic++) {
            if (ic==1) sspv = spectreV2;
            for(int jji=1;jji<sspv.Size()-1;jji++) 
              if (sspv(jji)>sspvmax[ic]) { sspvmax[ic]=sspv(jji); sspvmaxidx[ic]=jji; }
            if (nbsig[ic] < 1) { moysig[ic]=sigsig[ic]=-1.;  }
            else {
              moysig[ic] /= (double)nbsig[ic];
              sigsig[ic] /= (double)nbsig[ic];
              sigsig[ic] -= (moysig[ic]*moysig[ic]);
              sigsig[ic] = sqrt(sigsig[ic]);
              cout << "===Voie " << ic << " Moy=" << moysig[ic] << " Sig=" << sigsig[ic]
                   << " MaxSpec Amp= " << sqrt(sspvmax[ic])/double(pq.DataSize()/2/2) 
                   << " Pos=" << sspvmaxidx[ic] << "  (NPts=" << nbsig[ic] << ")" << endl;
            }
          }
          cout << "=== Voie1x2 MaxSpec Amp= " << sqrt(sqrt(sspvmax[2])/double(pq.DataSize()/2/2))
               << " Pos=" << sspvmaxidx[2] << endl;
        }  // fin if (fghist_)

        if (fgtimfreq) {
          timfreqV1 /= (r_4)nzm;
          timfreqV2 /= (r_4)nzm;
          po << PPFNameTag(tagtf1) << timfreqV1; 
          po << PPFNameTag(tagtf2) << timfreqV2; 
        }
        }
        spectreV1 = (r_4)(0.);
        spectreV2 = (r_4)(0.);
        visiV12 = complex<r_4>(0., 0.);
        if (fghist_) {
          ph1->Zero();
          ph2->Zero();
          moysig[0]=moysig[1]=0.;
          sigsig[0]=sigsig[1]=0.;
          nbsig[0]=nbsig[1]=0;
        }
        if (fgtimfreq) {
          timfreqV1 = (r_4)(0.);
          timfreqV2 = (r_4)(0.);
        }
        nzm = 0;  ifile++; 
//        ts.SetNow();
//        filog << ts << " :  proc file  " << fname << endl;                   
        cout << " BRProcA2C::run() created file  " << fname << card2name_(card_) << endl;
      }   
      
      memgr.FreeMemZone(mid, MemZS_ProcA);
    }  // Fin de boucle sur les zones a traiter 
  cout << " ------------  BRProcA2C::run() END " << card2name_(card_) 
       << " ------------ " << endl;
/*---- DELETE
  {
  nt.Info()["FirstTT"]=firsttt;     
  cout << nt;
  sprintf(fname,"%s_nt.ppf",path_.c_str());
  POutPersist po(fname);
  po << PPFNameTag("ntv12") << nt; 
  cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl;
  }
---- */ 
  if (nmoyspec>0) {  // Calcul des moyennes et sigmas des spectres 
    r_4 fnms = nmoyspec;
    moyspecV1 /= fnms;
    moyspecV2 /= fnms;
    sigspecV1 /= fnms;
    sigspecV2 /= fnms;
    sigspecV1 -= (moyspecV1 && moyspecV1);
    sigspecV2 -= (moyspecV2 && moyspecV2);
    sigspecV1 = Sqrt(sigspecV1);
    sigspecV2 = Sqrt(sigspecV2);
    TVector<r_4> rsbV1, rsbV2;   // Rapport signal/bruit
    moyspecV1.DivElt(sigspecV1, rsbV1, false, true);
    moyspecV2.DivElt(sigspecV2, rsbV2, false, true);
    sprintf(fname,"%s_ms.ppf",path_.c_str());
    POutPersist po(fname);
    po << PPFNameTag("moyspecV1") << moyspecV1; 
    po << PPFNameTag("moyspecV2") << moyspecV2; 
    po << PPFNameTag("sigspecV1") << sigspecV1; 
    po << PPFNameTag("sigspecV2") << sigspecV2; 
    po << PPFNameTag("rsbV1") << rsbV1; 
    po << PPFNameTag("rsbV2") << rsbV2; 
    cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl;
  }

  if (fghist_) {
    delete ph1; 
    delete ph2;
  }
  ts.SetNow();
  tm.SplitQ();
  cout << "  TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= " 
       << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" 
       << " ProcDataOut=" <<  totnbytesout/(1024*1024) << " MB" << endl;    
  cout << pcheck;
  cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl; 
  tm.Print();
  cout << " ---------------------------------------------------------- " << endl;
    
  }
  catch (PException& exc) {
    cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl;
    setRC(3);   
    return; 
  }
  catch(...) {
    cout << " BRProcA2C::run()/catched unknown ... exception " << endl;
    setRC(4);   
    return; 
  }
  setRC(0);
  return;
}   


//---------------------------------------------------------------------
// Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C)
//         !!!! OBSOLETE  !!!!   
//---------------------------------------------------------------------
/*!
   \class BRProcB4C
   \ingroup TAcq

   \brief Deprecated processing class. 
*/

/* --Methode-- */
BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path,
                     bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl)
  :  memgr1(mem1), memgr2(mem2) 
{
  fgraw_ = fgraw;
  nmax_ = nmax; 
  nmean_ = nmean;
  if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl;
  else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl;
  stop_ = false;        
  path_ = path; 
  fgnotrl_ = fgnotrl;
}

/* --Methode-- */
void BRProcB4C::Stop()
{
 stop_=true;
 // cout <<" BRProcB4C::Stop ... > STOP " << endl;
}


/* --Methode-- */
void BRProcB4C::run()
{
  setRC(1);     
  try {
    Timer tm("BRProcB4C", false);
    TimeStamp ts; 
    BRPaqChecker pcheck1(!fgnotrl_);  // Verification/comptage des paquets 
    BRPaqChecker pcheck2(!fgnotrl_);  // Verification/comptage des paquets 

    size_t totnbytesout = 0;
    size_t totnbytesproc = 0;

    cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ 
         << " NMean=" << nmean_ << endl;        
    cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl;

    uint_4 paqsz = memgr1.PaqSize();
    uint_4 procpaqsz = memgr1.ProcPaqSize();
    if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
      cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
           << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1=" 
           << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl;
      setRC(9); 
      return;   
    }

    TVector< complex<r_4> > cfour;  // composant TF
    BRPaquet pq(NULL, NULL, paqsz);     
/*   
    TVector<r_4> vx(pq.DataSize()/2);
    vx = (r_4)(0.);
    FFTPackServer ffts;
    ffts.FFTForward(vx, cfour);
    
    TVector< complex<r_4> > visiV13( cfour.Size() );
    TVector< complex<r_4> > visiV14( cfour.Size() );
    TVector< complex<r_4> > visiV23( cfour.Size() );
    TVector< complex<r_4> > visiV24( cfour.Size() );
*/
    int szfour = pq.DataSize()/2/2+1;
//    int szfour = (paqsz-40)/2+1;
    TVector< complex<r_4> > visiV13( szfour  );
    TVector< complex<r_4> > visiV14( szfour );
    TVector< complex<r_4> > visiV23( szfour );
    TVector< complex<r_4> > visiV24( szfour );
    // cout << " *DBG*AAAAA ---- Vectors OK" << endl; 
    cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz 
         << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
         << " *8="  << szfour*8 << endl;

    DataTable dt;
    dt.AddLongColumn("fc1");
    dt.AddLongColumn("tt1");
    dt.AddLongColumn("fc2");
    dt.AddLongColumn("tt2");
    DataTableRow dtr = dt.EmptyRow();

    uint_4 nzm = 0;
    uint_4 totnoksfc = 0;
    uint_4 totnokpaq = 0;
    uint_4 totnpaq = 0;
    uint_4 ifile = 0;

    uint_4 curfc=0;
    uint_8 curtt=0;
    uint_4 curfc2=0;
    uint_8 curtt2=0;
    uint_8 firsttt=0;
    uint_8 firsttt2=0;
    bool fgfirst=true;
    for (uint_4 kmz=0; kmz<nmax_; kmz++) {
      uint_4 noksfc = 0;
      uint_4 nokpaq = 0;
      if (stop_) break;
       // cout << " *DBG*BBBBB" << kmz << endl; 

      int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
      Byte* buff1 = memgr1.GetMemZone(mid1);
      if (buff1 == NULL) {
         cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
             break;             
      }
      Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
      if (procbuff1 == NULL) {
            cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
        break;  
      }
      int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
      Byte* buff2 = memgr2.GetMemZone(mid2);
      if (buff1 == NULL) {
         cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
             break;             
      }
      Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
      if (procbuff2 == NULL) {
            cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
        break;  
      }
      uint_4 i1,i2;
      i1=i2=0;
//      cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
      while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) {
        BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz); 
        BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz); 
        totnpaq++;
//        cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
//<<","<<paq2.FrameCounter()<<endl;
        // on ne traite que les paquets OK
        if (!pcheck1.Check(paq1)) { i1++; continue; } 
        if (!pcheck2.Check(paq2)) { i2++; continue; }
        nokpaq++;
        if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; }  
        if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; }  
//        cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
// <<","<<paq2.FrameCounter()<<endl;

        if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
          cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
          break;
        }
        // Les deux framecounters sont identiques ...
        noksfc++;
        curfc=paq1.FrameCounter();
        curfc2=paq2.FrameCounter();
        if (fgfirst) { 
          firsttt=paq1.TimeTag();  firsttt2=paq2.TimeTag(); 
          cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT=" 
               << firsttt<<" , "<<firsttt2 <<endl;
          fgfirst=false; 
        } 
        curtt=paq1.TimeTag()-firsttt;
        curtt2=paq2.TimeTag()-firsttt2;
        dtr[0]=curfc;  dtr[1]=curtt;
        dtr[2]=curfc2;  dtr[3]=curtt2;
        dt.AddRow(dtr);

        complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz);
        complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2);
        complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz);
        complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2);
        for(sa_size_t j=0; j<visiV13.Size(); j++) { 
          visiV13(j)+=zp1[j]*conj(zp3[j]);
          visiV14(j)+=zp1[j]*conj(zp4[j]);
          visiV23(j)+=zp2[j]*conj(zp3[j]);
          visiV24(j)+=zp2[j]*conj(zp4[j]);
        }
        if (nzm==0) {
          visiV13.Info()["StartFC"] = curfc;
          visiV14.Info()["StartFC"] = curfc;
          visiV23.Info()["StartFC"] = curfc;
          visiV24.Info()["StartFC"] = curfc;
          visiV13.Info()["StartTT"] = curtt;
          visiV14.Info()["StartTT"] = curtt;
          visiV23.Info()["StartTT"] = curtt;
          visiV24.Info()["StartTT"] = curtt;
        }
        nzm++;  i1++; i2++;
        totnbytesproc += 2*paq1.DataSize();
      } // Fin de boucle sur les paquets d'une zone
      memgr1.FreeMemZone(mid1, MemZS_ProcB);
      memgr2.FreeMemZone(mid2, MemZS_ProcB);

      if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
        visiV13 /= complex<r_4>((r_4)nzm, 0.);
        visiV14 /= complex<r_4>((r_4)nzm, 0.);
        visiV23 /= complex<r_4>((r_4)nzm, 0.);
        visiV24 /= complex<r_4>((r_4)nzm, 0.);
        visiV13.Info()["NPaqMoy"] = nzm;
        visiV14.Info()["NPaqMoy"] = nzm;
        visiV23.Info()["NPaqMoy"] = nzm;
        visiV24.Info()["NPaqMoy"] = nzm;
        visiV13.Info()["EndFC"] = curfc;
        visiV14.Info()["EndFC"] = curfc;
        visiV23.Info()["EndFC"] = curfc;
        visiV24.Info()["EndFC"] = curfc;
        visiV13.Info()["EndTT"] = curtt;
        visiV14.Info()["EndTT"] = curtt;
        visiV23.Info()["EndTT"] = curtt;
        visiV24.Info()["EndTT"] = curtt;
        char fname[512];
        {
        sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
        POutPersist po(fname);
        po << PPFNameTag("visiV13") << visiV13; 
        po << PPFNameTag("visiV14") << visiV14; 
        po << PPFNameTag("visiV23") << visiV23; 
        po << PPFNameTag("visiV24") << visiV24; 
        }
        visiV13 = complex<r_4>(0., 0.);
        visiV14 = complex<r_4>(0., 0.);
        visiV23 = complex<r_4>(0., 0.);
        visiV24 = complex<r_4>(0., 0.);
        nzm = 0;  ifile++; 
//        ts.SetNow();
//        filog << ts << " :  proc file  " << fname << endl;                   
        cout << " BRProcB4C::run() created file  " << fname << endl;
      }   
      double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
      cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc 
           << " (" << okfrac << " %)" << endl;
      totnokpaq += nokpaq;
      totnoksfc += noksfc;
    }  // Fin de boucle sur les zones a traiter 
    cout << " ------------------  BRProcB4C::run() END ----------------- " << endl;
    {     
    dt.Info()["FirstTT1"]=firsttt;     
    dt.Info()["FirstTT2"]=firsttt2;     
    cout << dt;
    char fname[512];
    sprintf(fname,"%s_fctt.ppf",path_.c_str());
    POutPersist po(fname);
    po << PPFNameTag("ttfc") << dt; 
    cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname  << endl;
    }     
    ts.SetNow();
    tm.SplitQ();
    cout << "  TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= " 
         << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl;
    double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.;
    cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC=" 
         << totnoksfc << " (" << totokfrac << " %)" << endl;
//  cout << pcheck1;
//  cout << pcheck2;
    cout << " BRProcB4C::run()/Timing: \n"; 
    tm.Print();
    cout << " ---------------------------------------------------------- " << endl;
}
  catch (PException& exc) {
    cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl;
    setRC(3);   
    return; 
  }
  catch(...) {
    cout << " BRProcB4C::run()/catched unknown ... exception " << endl;
    setRC(4);   
    return; 
  }
  setRC(0);
  return;
}