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

//----------------------------------------------------------------
// Projet BAORadio - (C) LAL/IRFU  2008-2010
// 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 de traitement simple - calcul de spectres moyennes / voie
//---------------------------------------------------------------------
/* --Methode-- */
BRMeanSpecCalculator::BRMeanSpecCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, 
                                           bool fgdatafft, bool fgsinglechan)
  : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean),
    fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan), 
    clnflg_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()),
    nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()) 
{
  setNameId("meanSpecCalc",1);
  BRPaquet paq(memgr_.PaqSize());
  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;
  nbpaq4mean_=0;
  totnbpaq_=0;

  SetVarianceLimits();

  ofsdtp_=NULL;
  dtp_=NULL;
  string dtfile="!"+outpath_+"/dtspec.fits";
  ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
  dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true);
  int nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
  char cnom[32];
  for(int i=0; i<nchan; i++) {
    sprintf(cnom,"var%d",i);
    dtp_->AddFloatColumn(cnom);    
  }
  for(int i=0; i<nchan; i++) {
    sprintf(cnom,"varnorm%d",i);
    dtp_->AddFloatColumn(cnom);    
  }
  xnt_=new double[nchan*2];
}

/* --Methode-- */
BRMeanSpecCalculator::~BRMeanSpecCalculator()
{
  if (nbpaq4mean_>1)  SaveSpectra();
  cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
  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_;
    delete xnt_;
  }
  cout << " ------------------------------------------------------------------------ " << endl;
}

/* --Methode-- */
void BRMeanSpecCalculator::ReadGainFitsFile(string filename)
{
  cout << " BRMeanSpecCalculator::ReadGainFitsFile() - reading file " << filename;
  FitsInOutFile fis(filename, FitsInOutFile::Fits_RO);
  fis >> sgain_;
  cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << endl;
}

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

/* --Methode-- */
int BRMeanSpecCalculator::Process()
{
  if (nbpaq4mean_==nmean_)  SaveSpectra();
  FlagBadPackets();
    
  if (fgdatafft_) {  // Donnees firmware FFT 
    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
      if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) 
      TwoByteComplex* zp=NULL;
      if (fgsinglechannel_) {
        zp=vpaq_[i].Data1C();
      }
      else {
        zp=vpaq_[i/2].Data1C();
        if (i%2==1)  zp=vpaq_[i/2].Data2C();
      }
      TVector< r_4 > spec = mspecmtx_.Row(i); 
      TVector< r_4 > sspec = sigspecmtx_.Row(i); 
      for(sa_size_t f=1; f<spec.Size(); f++) {
        r_4 m2zf=zp[f].module2F();;
        spec(f) += m2zf;
        sspec(f) += m2zf*m2zf;
      }
    }
  }
  else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator
    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
      if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) 
      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) ;
      }
      TVector< r_4 > spec = mspecmtx_.Row(i); 
      TVector< r_4 > sspec = sigspecmtx_.Row(i); 
      for(sa_size_t f=1; f<spec.Size(); f++)  {
        r_4 m2zf=Zmod2(zp[f]);
        spec(f) += m2zf;
        sspec(f) += (m2zf*m2zf);        
      }     
    }
  }
  nbpaq4mean_++;   totnbpaq_++;
  return 0;
}

/* --Methode-- */
void BRMeanSpecCalculator::FlagBadPackets()
{
  if (fgdatafft_) {  // Donnees firmware FFT 
    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); 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();
      }
      double variance=0.;
      double varnorm=0.;
      for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) {
        double modsq=(double)zp[f].module2F();
        variance+=modsq;
        varnorm+=modsq/sgain_(i,f);
      }
      xnt_[i]=variance;
      xnt_[i+mspecmtx_.NRows()]=varnorm;
      clnflg_[i]=0;
      if (varnorm<varmin_) clnflg_[i]=1;
      else if (varnorm>varmax_) clnflg_[i]=2;
    }
  }
  else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator
    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); 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) ;
      }
      double variance=0.;
      double varnorm=0.;
      for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) {
        double modsq=(double)Zmod2(zp[f]);
        variance+=modsq;
        varnorm+=modsq/sgain_(i,f);
      }
      xnt_[i]=variance;
      xnt_[i+mspecmtx_.NRows()]=varnorm;
      clnflg_[i]=0;
      if (varnorm<varmin_) clnflg_[i]=1;
      else if (varnorm>varmax_) clnflg_[i]=2;
    }
  }
  dtp_->AddRow(xnt_);
  return;
}

/* --Methode-- */
void BRMeanSpecCalculator::SaveSpectra()
{
  mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_;
  mspecmtx_ /= (double)nbpaq4mean_;
  sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_;
  sigspecmtx_ /= (double)nbpaq4mean_;
  sigspecmtx_ -= (mspecmtx_ && mspecmtx_);  // Mean(X^2) - [ Mean(X) ]^2
  char nfile[64];
  string flnm;
  /*
  {
  sprintf(nfile,"mspecmtx%d.ppf",numfile_);
  flnm=outpath_+nfile;
  POutPersist po(flnm);
  po << mspecmtx_;
  }
  {
  sprintf(nfile,"sigspecmtx%d.ppf",numfile_);
  flnm=outpath_+nfile;
  POutPersist po(flnm);
  po << sigspecmtx_;
  }
  */
  {
  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::SaveSpectra() NPaqProc="
       << totnbpaq_ << "  -> Mean/Sig spectra Matrix in " << flnm << " /sigspec...ppf" << endl;
  numfile_++;  nbpaq4mean_=0;
  mspecmtx_ = (r_4)(0.);
  sigspecmtx_ = (r_4)(0.);
  return;
}

//---------------------------------------------------------------------
// Classe de thread de calcul de FFT sur donnees RAW 
//---------------------------------------------------------------------
/* --Methode-- */
BRFFTCalculator::BRFFTCalculator(RAcqMemZoneMgr& memgr, bool fgsinglechannel)
  : BRBaseProcessor(memgr), fgsinglechannel_(fgsinglechannel), totnbfftpaq_(0)
{
  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)
//-------------------------------------------------------------------------
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);
  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  !!!!   
//---------------------------------------------------------------

// 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) += 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) += Zmod2(cfour1(j));
          }
        }
        for(sa_size_t j=0; j<spectreV2.Size(); j++) 
          spectreV2(j) += Zmod2(cfour2(j));  // 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) += 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 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j)); 
          complex<r_4> zvis =  cfour1(j)*conj(cfour2(j));
          ms12 += 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 = 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  !!!!   
//---------------------------------------------------------------------

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