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

#include "racquproc.h"

#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 "fftwserver.h"

#include "FFTW/fftw3.h"


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

//---------------------------------------------------------------------
// Classe de traitement - calcul de visibilite pour n fibres 
//---------------------------------------------------------------------

/* --Methode-- */
BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, 
                                               uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
  : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean), numfreq1_(freq1), numfreq2_(freq2), nbinfreq_(nbfreq),
    vpdata_(2*memgr.NbFibres())
    // , dtfos_(outpath+"visdt.fits", Fits_Create), visdt_(dtfos_, 1024, true);

{  
  BRPaquet paq(memgr_.PaqSize());

  vismtx_.SetSize((2*memgr_.NbFibres()+1)*memgr_.NbFibres(), paq.DataSize()/4); 
  chanum_.SetSize(vismtx_.NRows());
  sa_size_t k=0;
  for(size_t i=0; i<2*memgr_.NbFibres(); i++)  vpdata_[i]=NULL;
  for(size_t i=0; i<2*memgr_.NbFibres(); i++) {
    for(size_t j=i; j<2*memgr_.NbFibres(); j++) {
      chanum_(k) = (i+1)*100+(j+1);  k++;
    }
  }
  //  visdt_.AddFloatColumn("mfc");
  visdt_.AddDoubleColumn("mfc");
  visdt_.AddDoubleColumn("mtt");
  visdt_.AddIntegerColumn("jfreq");
  visdt_.AddIntegerColumn("numch");
  visdt_.AddFloatColumn("vre");
  visdt_.AddFloatColumn("vim");

  if (nmean_ < 1)  nmean_=memgr_.NbPaquets();
  if (nmean_ < 1)  nmean_=1;

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

  jf1_=numfreq1_; jf2_=numfreq2_;
  if ((jf1_<1)||(jf1_>=vismtx_.NCols()))  jf1_=1;
  if ((jf2_<1)||(jf2_>=vismtx_.NCols())||(jf2_<jf1_))  jf2_=vismtx_.NCols()-1;
  if (nbinfreq_<1) nbinfreq_=1;
  djf_=(jf2_-jf1_)/nbinfreq_;
  if (djf_<1) djf_=0;
  xnt_ = new double[vismtx_.NRows()+20];
  for(int kk=0; kk<vismtx_.NRows(); kk++)  xnt_[kk]=0.;

  cout << " BRVisibilityCalculator/Info NMean= " << nmean_ << " JF1=" << jf1_ 
       << " JF2=" << jf2_ << " DJF=" << djf_ << endl;

  fgcktt_=false;
}

/* --Methode-- */
BRVisibilityCalculator::~BRVisibilityCalculator()
{
  cout << " BRVisibilityCalculator - Visibility Datatable : " << endl;
  cout << visdt_;
  POutPersist po(outpath_+"visdt.ppf");
  po << visdt_;
  POutPersist poc(outpath_+"chanum.ppf");
  poc << chanum_;

  if (fgcktt_) {
    cout << " BRVisibilityCalculator -  Check TimeTag done: " << endl;
    for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  {
      cout << " BRTTCheck-Fiber[" << fib << "] NBadTT=" << vbadtt_[fib] 
           << "  NotSameTT=" << vnsamett_[fib] << endl; 
    }
    POutPersist pott(outpath_+"ttfcmtx.ppf");
    pott << PPFNameTag("FC") << fcmtx_; 
    pott << PPFNameTag("TT") << ttmtx_;     
  }
  delete[] xnt_;
}

/* --Methode-- */
int BRVisibilityCalculator::ActivateTimeTagCheck(uint_8 maxnpaq)
{
  mindeltatt_=memgr_.PaqSize()/2;
  if (mindeltatt_<1) mindeltatt_=1;
  fcmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
  ttmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
  vlasttt_.resize(memgr_.NbFibres(), 0);
  vbadtt_.resize(memgr_.NbFibres(), 0);
  vnsamett_.resize(memgr_.NbFibres(), 0);
  fgcktt_=true;
  cout << " BRVisibilityCalculator::ActivateTimeTagCheck() - TT/Fc matrix NCols=" << maxnpaq 
       << " MinDeltaTT=" << mindeltatt_ << endl;

  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::Process()
{
    
  for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
    vpdata_[2*fib] = vpaq_[fib].Data1C();
    vpdata_[2*fib+1] = vpaq_[fib].Data2C();
  }
  
  if (totnbpaq_%nmean_ == 0) {
    if (totnbpaq_ > 0) {
      moyfc_/=nmean_;
      moytt_/=nmean_;
      vismtx_.Info()["MeanFC"] = moyfc_;
      vismtx_.Info()["MeanTT"] = moytt_;
      vismtx_.Info()["NPAQSUM"] = nmean_;

      // ATTENTION : Matrice visibilites non moyennee 
      char nfile[32];
      sprintf(nfile,"vismtx%d.ppf",numfile_);
      string flnm=outpath_+nfile;
      POutPersist po(flnm);
      po << vismtx_;
      cout << numfile_ << "-BRVisibilityCalculator::Process() NPaqProc="
           << totnbpaq_ << "  -> Visibility Matrix in " << flnm << endl;
      FillVisibTable(moyfc_, moytt_);
      numfile_++;
    }
    vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
    moyfc_=moytt_=0.;
  }
  
  sa_size_t k=0;
  for(size_t i=0; i<vpdata_.size(); i++) {
    for(size_t j=i; j<vpdata_.size(); j++) {
      TVector< complex<r_4> > vis = vismtx_.Row(k);   k++;
      for(sa_size_t f=1; f<vis.Size(); f++) {
        vis(f) += complex<r_4>((r_4)vpdata_[i][f].realB(), (r_4)vpdata_[i][f].imagB()) * 
          complex<r_4>((r_4)vpdata_[j][f].realB(), -(r_4)vpdata_[j][f].imagB());
      }
    }
  }
  
  moyfc_ += curfc_[0];
  moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
  if (fgcktt_)  CheckTimeTag();
  totnbpaq_++;
  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
{
  xnt_[0]=fcm;  xnt_[1]=ttm/1.25e8;

  if (djf_<2) {
    for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
      for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
        xnt_[2]=jf;
        xnt_[3]=chanum_(rv);
        xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
        xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
        visdt_.AddRow(xnt_);
      }
    }
  }
  else {
    for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
      for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
        r_4 moyreal=0.;
        r_4 moyimag=0.;
        sa_size_t jjfmx=jf+djf_;
        if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
        for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
          moyreal+=vismtx_(rv,jjf).real();
          moyimag+=vismtx_(rv,jjf).imag();
        }
        xnt_[2]=jf+djf_/2;
        xnt_[3]=chanum_(rv);
        xnt_[4]=moyreal/(r_4)(nmean_*djf_);
        xnt_[5]=moyimag/(r_4)(nmean_*djf_);
        visdt_.AddRow(xnt_);
      }
    }
  }
  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::CheckTimeTag()
{
  if (totnbpaq_==0) {
    for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  {
      vlasttt_[fib]=ttfirst_[fib];
      if (ttmtx_.NCols()>0) {
        fcmtx_(fib,totnbpaq_) = curfc_[fib]; 
        ttmtx_(fib,totnbpaq_) = vlasttt_[fib]; 
      }
    }
    return 0;
  }
  for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
    uint_8 ld = vpaq_[fib].TimeTag()-vlasttt_[fib];
    uint_8 fd = vpaq_[fib].TimeTag()-ttfirst_[fib]-vpaq_[0].TimeTag()+ttfirst_[0];
    /*    if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
    else {
      if (fd!=0)  vnsamett_[fib]++;
    }
    */
    if  (ld < mindeltatt_) vbadtt_[fib]++;
    if (fd != 0)  vnsamett_[fib]++;
    vlasttt_[fib]=vpaq_[fib].TimeTag();
    if (totnbpaq_<ttmtx_.NCols()) {
      fcmtx_(fib,totnbpaq_) = curfc_[fib]; 
      ttmtx_(fib,totnbpaq_) = vlasttt_[fib]; 
    }
  }
  return 0;
}

//---------------------------------------------------------------
// Classe thread de traitement donnees ADC avec 2 voies par frame
//---------------------------------------------------------------

// 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 r_4 Zmod2(complex<r_4> z) 
{ return (z.real()*z.real()+z.imag()*z.imag()); }

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

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