source: Sophya/trunk/AddOn/TAcq/corrcrtadc.cc@ 4020

// Utilisation de SOPHYA ...
#include "sopnamsp.h"
#include "machdefs.h"

/* --------------------------------------------------------------- 
   Projet BAORadio-21cm intensity mapping - (C) LAL/IRFU  2008-2011

   CRT (Pittsburgh-CMU) DUMP ADC read/processing program to compute 
   spectra and correlations du CRT
   R. Ansari, C. Magneville   -  LAL/Irfu - Juin 2011
   --------------------------------------------------------------- */

// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <string>

// SOPHYA include files 
#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h" 
#include "datatable.h" 
#include "histinit.h" 
#include "matharr.h" 
#include "timestamp.h"
#include "fftwserver.h"
#include "fftpserver.h"
#include "utilarr.h"
#include "histats.h"

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"


//--------------------------- Functions in this file  ------------------- 
int Usage(bool fgshort=true);
int DecodeProc(int narg, char* arg[]);
int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname, bool fgcombiq, 
                    int imin, int imax, int istep,  int jf1, int jf2, int nfreq);
int ComputeCorrel(TMatrix< complex<r_4> >& cfour,  bool fgcombiq, 
                  TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel);
int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep);
int CombineIQPairs(TMatrix< complex<r_4> >& cfour);
//------------------------------------------------------------------------------------------------------------

/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- corrcrtadc.cc : Read CRT-ADC dump files and process to produce correlation matrices \n"
       << " Usage: corrcrtadc -corr/-corciq/-hval InOutPath OutFile [AdcIds] [Imin,Imax,step]\n" << endl;
  if (fgshort) {
    cout << " corrcrtadc -h for detailed instructions" << endl;
    return 1;
  }
  cout << " -corr/-corciq/-hval : Compute correlations or Fill histograms for ADC sample values \n"
       << "    -corciq : combine IQ pairs to separate side-bands ... \n" 
       << " InOutPath : Input/Output directory name \n" 
       << " ADCIds : ADC Id's to be processed (Example: 1,2,3,4, default=0,1)  \n" 
       << " OutFile : Output PPF file  name \n" 
       << " Imin,Imax,IStep: Input ADC dump files sequence number \n"
       << "    FileNames=InOutPath/adcdumpNJFII.fits Imin<=II<=Imax II+=IStep , J=ADCIds \n" << endl;
  //       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple" << endl;
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<4) return Usage(true);
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
    DecodeProc(narg-1, arg+1);
    resu.Update();
    cout << resu;
  }
  catch (PException& exc) {
    cerr << " corrcrtadc.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n corrcrtadc.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " corrcrtadc.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> corrcrtadc.cc ------- END ----------- RC=" << rc << endl;
  return rc;

}


//-------------------------------------------------------------------- 
//       Traitement fichiers produits par vismfib  (V Nov09)
/* --Fonction-- */
int DecodeProc(int narg, char* arg[])
{
  // Decodage des arguments et traitement 
  string popt = arg[0];
  bool fillhis=false;
  bool fgcomiq=false;
  if (popt=="-hval")  fillhis=true;
  else if (popt=="-corciq")  fgcomiq=true;
  string inoutpath = arg[1];  
  string outname = arg[2];
  int aids[8]={0,1,-1,-1,-1,-1,-1,-1};
  if (narg>3) sscanf(arg[3],"%d,%d,%d,%d,%d,%d,%d,%d",aids,aids+1,aids+2,aids+3,aids+4,aids+5,aids+6,aids+7);
  vector<int> adcids;
  for(int ii=0;ii<8;ii++) 
    if((aids[ii]>=0)&&(aids[ii]<16))  adcids.push_back(aids[ii]);
  int imin=0;
  int imax=0;
  int istep=1;
  if (narg>4) sscanf(arg[4],"%d,%d,%d",&imin,&imax,&istep);
  int jf1=0;
  int jf2=0;
  int nfreq=0;
  if (narg>5)
    sscanf(arg[5],"%d,%d,%d",&jf1,&jf2,&nfreq);

  cout << " ----- corrcrtadc/DecodeProc - Start " << ((fillhis)? " Fill ADC-val histogram" : " Compute Correlation") << endl;
  cout <<  " InOutPath= " << inoutpath << " ADCIds: " ;
  for(size_t ii=0; ii<adcids.size(); ii++) cout << adcids[ii] << " , ";
  cout << " (NbADC=" << adcids.size() << ") OutFileName= " << outname 
       << " IMin,Max,Step="    << imin << "," << imax << "," << istep << endl;
  //       << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
  int rc=0;
  if (fillhis) rc = ADCFiles2Histo(inoutpath, adcids, outname, imin, imax, istep);
  else rc=ProcessADCFiles(inoutpath, adcids, outname, fgcomiq, imin, imax, istep, jf1, jf2, nfreq);
  
  return rc;
}

// Taille des fichiers 4 voies = 1024*1024*128 ===> 32 x 1024 x 1024 samples / voie
// On lit par paquet de 4096 = 4 x 1024 / voie --> ADCRDBLKSZ=16x1024=16384
#define ADCFLEN 4096 
#define ADCRDBLKSZ 16384
#define ADCNREAD 8192 

static int PRTMODULO = 1000;

// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
/* --Fonction-- */
int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname,  bool fgcombiq, int imin, int imax, int istep, int jf1, int jf2, int nfreq) 
{
  Timer tm("ProcessADCFiles");  

#define NFREQUSED 2040
#define NFREQUSESTART 5
#define NFREBIN 4 

  sa_size_t NFREQREBIN = NFREQUSED/NFREBIN;
  sa_size_t NBADC = adcids.size();

  sa_size_t  NCHAN = 4*NBADC;
  TMatrix< complex<r_4> > cfour(NCHAN, NFREQUSED);
  sa_size_t  NCORREL = NCHAN*(NCHAN+1)/2;
  TMatrix< complex<r_8> > correl(NCORREL,NFREQUSED);
  TMatrix< complex<r_8> > autocorrel(NCHAN,NFREQUSED);
  double nsumcor=0.;

  vector<FILE*> fip; 
  vector<int> adcchans;
  for(size_t ka=0; ka<NBADC; ka++)  { 
    fip.push_back(NULL);
    for(size_t ica=0; ica<4; ica++) 
      adcchans.push_back(adcids[ka]*4+ica);
  }

  FFTPackServer ffts(false);
  
  TVector<r_4> sigadc(ADCFLEN);
  TVector< complex<r_4> > foursig;
  
  char fname[512];

  signed char rdbuff[ADCRDBLKSZ]; 
  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    for(int ka=0; ka<NBADC; ka++) {
      sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
      fip[ka]=fopen(fname,"rb");
      if (fip[ka]==NULL) {
        cout << "ProcessADCFiles/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
        throw IOExc(" ADC dump file open error!");
      }
      else  
        cout << "ProcessADCFiles/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
    }

    for(int m=0; m<ADCNREAD; m++) {
      sa_size_t irc=0;
      for(int ka=0; ka<NBADC; ka++) {
        fread(rdbuff,1,ADCRDBLKSZ,fip[ka]);
        for(int jac=0; jac<4; jac++) {
          for(sa_size_t ls=0; ls<ADCFLEN; ls++) sigadc(ls)=(float)rdbuff[ls*4+jac];
          ffts.FFTForward(sigadc, foursig);
          cfour.Row(irc)=foursig(Range(NFREQUSESTART,NFREQUSESTART+NFREQUSED-1)).Transpose();  irc++;
        }
      }
      ComputeCorrel(cfour, fgcombiq, correl, autocorrel);
      nsumcor++;
      if (m%PRTMODULO==0)  cout << " ProcessADCFiles/Info Done read+FFT+correl m=" << m << " / Max=" << ADCNREAD << endl;
    }
    for(int ka=0; ka<NBADC; ka++)  fclose(fip[ka]);
  }


  correl /= nsumcor;
  autocorrel /= nsumcor;

  TMatrix< complex<r_8> > correbin(NCORREL,NFREQREBIN);
  sa_size_t icc=0;
  for(sa_size_t i=0; i<correbin.NCols(); i++) {
    for(sa_size_t j=0; j<NFREBIN; j++) { 
      correbin.Column(i) += correl.Column(icc); 
      icc++;
    }
  }

  sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
  cout << "ProcessADCFiles: Opening file " << fname << " for writing correlation matrix" << endl;  
  POutPersist po(fname);
  po << PPFNameTag("correlmtx") << correl; 
  po << PPFNameTag("acorrmtx") << autocorrel; 
  po << PPFNameTag("rebincormtx") << correbin; 

  cout << "ProcessADCFiles:  finished FFT + correlation computing " << endl;

  cout << " --------------- Channel/Correlation number association --------------- " << endl;
  cout << " CorreNumber- axb  (MtxRow) :  (Ca, Cb)  -- ADCChan(c,d) -- I/Q(e,f)" << endl; 
  cout << " ----------------------------------------------------------------------- " << endl;
  sa_size_t jcr=0;
  char IQid[2];
  char EWid[2];
  int EWnum[2];
  int Bid[2];
  int Inpid[2];

  for(sa_size_t j1=0; j1<cfour.NRows(); j1++) {
    for(sa_size_t j2=j1; j2<cfour.NRows(); j2++) {
      cout <<  jcr << " : (" << j1 << "," << j2 << ") -> ADCChans (" << adcchans[j1] << "," << adcchans[j2] << ") -- (";
      Bid[0]=adcchans[j1]/4;
      Bid[1]=adcchans[j2]/4;
      Inpid[0]=adcchans[j1]%4;
      Inpid[1]=adcchans[j2]%4;

      if ((adcchans[j1]/4)%2==0) IQid[0]='I'; else IQid[0]='Q'; 
      if ((adcchans[j2]/4)%2==0) IQid[1]='I'; else IQid[1]='Q'; 
      if (adcchans[j1]%2==0) EWid[0]='E'; else EWid[0]='W'; 
      if (adcchans[j2]%2==0) EWid[1]='E'; else EWid[1]='W'; 
      cout << IQid[0] << '-' << EWid[0] << (Bid[0]/2)*2+Inpid[0]/2 << "," 
           << IQid[1] << '-' << EWid[1] << (Bid[1]/2)*2+Inpid[1]/2 << ")" << endl;
      jcr++;
    }
  }
  cout << " ----------------------------------------------------------------------- " << endl;

  return 0;
}

/* --Fonction-- */
int ComputeCorrel(TMatrix< complex<r_4> >& cfour,  bool fgcombiq, 
                  TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel)
{
  if (fgcombiq)  CombineIQPairs(cfour);
  sa_size_t jcr=0;
  for(sa_size_t j1=0; j1<cfour.NRows(); j1++) {
    for(sa_size_t j2=j1; j2<cfour.NRows(); j2++) {
      for(sa_size_t i=0; i<cfour.NCols(); i++) 
        correl(jcr,i) += complex<r_8> ( cfour(j1,i)*conj(cfour(j2,i)) ); 
      if (j1==j2) autocorrel.Row(j1)=correl.Row(jcr);
      jcr++;
    }
  }
  return 0;
}

/* --Fonction-- */
int CombineIQPairs(TMatrix< complex<r_4> >& cfour)
{
  // We assume the following channel numbering for the I-Q pairs
  // Row pairs : (0,4) (1,5) (2,6) (3,7) (8,12) (9,13) (10,14) (11,15) ... are assumed to be 
  // IQ pairs in the matrix cfour 
  // We compute the I rows by the left side band, and the Q rows by the right side band 

  // Left = a I + b Q
  // Right = c I + d Q 
  complex<r_4> za(1.,0.);
  complex<r_4> zb(0.,1.);
  complex<r_4> zc(1.,0.);
  complex<r_4> zd(0.,-1.);

  int nboards = cfour.NRows()/4;    
  for(int jb=0; jb<nboards-1; jb++)  {
    for(int c=0; c<4; c++) {
      sa_size_t ii = jb*4+c;  // sa_size_t are long integers 
      sa_size_t iq = ii+4;
      //  Combine the two I,Q fourier coefficient vectors to form the left and right sid bands 
      TVector< complex<r_4> > leftband  = za*cfour.Row(ii)+zb*cfour.Row(iq); 
      TVector< complex<r_4> > rightband = zc*cfour.Row(ii)+zd*cfour.Row(iq); 
      // replace the two I,Q vectors by left and right computed side bands 
      cfour.Row(ii)=leftband;
      cfour.Row(iq)=rightband;
    }
  }
  return 0;
}

/* --Fonction-- */
int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep)
{
  Timer tm("ADCFiles2Histo");  

  sa_size_t NBADC = adcids.size();

  sa_size_t  NCHAN = 4*NBADC;
  vector<Histo *> vhist;
  for(int ka=0; ka<NCHAN; ka++) 
    vhist.push_back(new Histo(-128.5,128.5,257) );

  vector<FILE*> fip; 
  for(size_t ka=0; ka<NBADC; ka++)  fip.push_back(NULL);
    
  char fname[512];

  signed char rdbuff[ADCRDBLKSZ]; 
  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    for(int ka=0; ka<NBADC; ka++) {
      sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
      fip[ka]=fopen(fname,"rb");
      if (fip[ka]==NULL) {
        cout << "ADCFiles2Histo/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
        throw IOExc(" ADC dump file open error!");
      }
      else  
        cout << "ADCFiles2Histo/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
    }

    for(int m=0; m<ADCNREAD; m++) {
      sa_size_t irc=0;
      for(int ka=0; ka<NBADC; ka++) {
        fread(rdbuff,1,ADCRDBLKSZ,fip[ka]);
        for(int jac=0; jac<4; jac++) {
          for(sa_size_t ls=0; ls<ADCFLEN; ls++) vhist[irc]->Add((r_8)rdbuff[ls*4+jac]);
          irc++;
        }
      }
      if (m%PRTMODULO==0)  cout << " ADCFiles2Histo/Info Done FillHist m=" << m << " / Max=" << ADCNREAD << endl;
    }
    for(int ka=0; ka<NBADC; ka++)  fclose(fip[ka]);
  }

  sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
  cout << "ADCFiles2Histo: Opening file " << fname << " for writing ADC value histograms " << endl;  
  POutPersist po(fname);
  char nametag[64];
  for(size_t ka=0; ka<NCHAN; ka++) {
    sprintf(nametag,"adcvalC%d",ka);
    po << PPFNameTag(nametag) << *(vhist[ka]); 
  }
  for(size_t ka=0; ka<NCHAN; ka++) delete vhist[ka];

  cout << "ADCFiles2Histo:  finished ADC value histogram fill " << endl;
  return 0; 
}