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

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio - (C) LAL/IRFU  2008-2010

   Programme de lecture des fichiers DUMP ADC du CRT (Pittsburgh) 
   pour calcul de spectres / correlations 
   R. Ansari, C. Magneville   -  LAL/Irfu - Juin 2011
   V : Mai 2009
   ---------------------------------------------------------- */

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

#include <iostream>
#include <string>

#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 sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"


//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
int DecodeProc(int narg, char* arg[]);
int ProcessADCFiles(string& inoutpath, string& outname, int imin, int imax, int istep,  int jf1, int jf2, int nfreq);
int ComputeCorrel(TMatrix< complex<r_4> >& cfour,  TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel);
//------------------------------------------------------------------------------------------------------------

/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- corrcrtadc.cc : Read CRT-ADC dump files to produce time-frequency\n"
       << "     correlation matrices " << endl;
  cout << " Usage: corrcrtadc InOutPath OutFile Imin,Imax,step [NumFreq1,NumFreq2,NBinFreq] \n" << endl;
  if (fgshort) {
    cout << " corrcrtadc -h for detailed instructions" << endl;
    return 1;
  }
  cout << " InOutPath : Input/Output directory name \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=0,1,2 \n" 
       << " 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 inoutpath = arg[0];  
  string outname = arg[1];
  int imin=0;
  int imax=0;
  int istep=1;
  sscanf(arg[2],"%d,%d,%d",&imin,&imax,&istep);
  int jf1=0;
  int jf2=0;
  int nfreq=0;
  if (narg>3)
    sscanf(arg[3],"%d,%d,%d",&jf1,&jf2,&nfreq);

  cout << " ----- corrcrtadc/DecodeProc - Start - InOutPath= " << inoutpath << " OutFileName=" << outname << endl;
  cout << " IMin,Max,Step="    << imin << "," << imax << "," << istep 
       << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
  int rc=ProcessADCFiles(inoutpath, outname, imin, imax, istep, jf1, jf2, nfreq);
  return rc;
}

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

  // 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 
#define NBADC 2
#define NFREQUSED 2000

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

  FILE* fip[NBADC]={NULL,NULL};
  FFTPackServer ffts(false);
  
  TVector<r_4> sigadc(ADCFLEN);
  TVector< complex<r_4> > foursig;
  
  char fname[512];

  signed char rdbuff[ADCRDBLKSZ]; 
  int prtmodulo=50;
  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(),ka,ifile);
      fip[ka]=fopen(fname,"rb");
      if (fip[ka]==NULL) cout << "ProcessADCFiles/ERROR opening file " << fname << endl;
      else  cout << "ProcessADCFiles/Info - file " << fname << " opened ..." << 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(400,1679)).Transpose();  irc++;
        }
      }
      ComputeCorrel(cfour, correl, autocorrel);
      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]);
  }


  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; 

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

/* --Fonction-- */
int ComputeCorrel(TMatrix< complex<r_4> >& cfour,  TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel)
{
  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;
}