source: Sophya/trunk/AddOn/TAcq/mcrd.cc@ 4005

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

/* ---------------------------------------------------------- 
   Programme de lecture multi canaux pour BAORadio  
   R. Ansari, C. Magneville
   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 "timestamp.h"
#include "fftpserver.h"
#include "fftwserver.h"

#include "FFTW/fftw3.h"

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


// include mini-fits lib , et structure paquet BAORadio
#include "minifits.h"
#include "brpaqu.h"
#include "brfitsrd.h"
#include "brproc.h"
#include "racqurw.h"


//--- Fonctions de ce fichier 
int DecodeMiniFitsHeader(string& filename, uint_4& paqsz, uint_4& npaq, bool fgnotrl=false) 
{
  MiniFITSFile mff(filename, MF_Read);
  cout << "DecodeMiniFitsHeader()... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
       << " NAxis2=" << mff.NAxis2() << endl;
  paqsz = mff.NAxis1();
  npaq = mff.NAxis2();
  if (fgnotrl)  paqsz += 16;
  return 0;
}

//--  Parametres globaux
static int NZones=4;
static int NPaqinZone=128;
static int NMean=1024;
static int NFreqSMap=0;  // binning en frequences cartes 2D temps-freq, 0-> NoMap 2D spectres   
static int NGenZ=100;
static int GPaqSz=16424;
static double LossRate=0.1;

static bool fgraw=true;  // true -> RAW data , false -> FFT data 
static bool fg4c=false;  // true -> 4 channels (2 fibers)
static bool fgrdfits=true;  // false -> Don't read fits files, generate paquets
static bool fgnotrl=false;   // true -> fichier fits SANS Trailer de frame (< mai 2009)
static bool fghist=false;   // true -> histo des valeurs des time sample

// ----
int Usage(bool fgshort=true);
// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
int Proc1FA(string& outname, string& inpath, int jf1, int jf2);
// Pour traitement (calcul FFT et visibilites (ProcA,ProcB) 2 fibre, 4 voies RAW)
int Proc2FAB(string& outname, string& path1, string& path2, int jf1, int jf2);

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

  TArrayInitiator  _inia;

  int rc = 0;
  try {
    string act = arg[1];
    fg4c=false;  // true -> 4 channels (2 fibers)
    fgrdfits=true;  // false -> Don't read fits files, generate paquets
    fgnotrl=false;   // true -> fichier fits SANS Trailer de frame (< mai 2009)
    fghist=false;   // true -> histo des valeurs des time sample
    fgraw=true;     // raw data 
    if (act.substr(0,2)=="-4") fg4c=true;
    if (act.length()>2) {
      for(size_t ks=2; ks<act.length(); ks++) {
        if (act[ks]=='f') fgraw=false;
        else if(act[ks]=='g') fgrdfits=false;
        else if(act[ks]=='n') fgnotrl=true; 
        else if(act[ks]=='h') fghist=true; 
      }
    }
    if (fg4c && (narg<6)) return Usage(true);

    string outname = arg[2];    
    string inpath2;
    string inpath = arg[3]; 
    int offa=4;
    if(fg4c) { 
      inpath2=arg[4];    offa=5; 
    } 
    int imin=0;
    int imax=0;
    sscanf(arg[offa],"%d,%d",&imin,&imax);
    if (narg>offa+1) sscanf(arg[offa+1],"%d,%d",&NZones,&NPaqinZone);   
    if (narg>offa+2) NMean=atoi(arg[offa+2]);
    if (narg>offa+3) NFreqSMap=atoi(arg[offa+3]);
    if (narg>offa+4) sscanf(arg[offa+4],"%d,%d,%lg",&GPaqSz,&NGenZ,&LossRate);  

    cout << " ---------- mcrd.cc Start - ACT= " << act << " ------------- " << endl;
    ResourceUsage resu;
    if (fg4c) 
      rc =  Proc2FAB(outname, inpath, inpath2, imin, imax); 
    else 
      rc =  Proc1FA(outname, inpath, imin, imax); 
    cout << resu ; 
  }
  catch (MiniFITSException& exc) {
    cerr << " mcrd.cc catched MiniFITSException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n mcrd.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " mcrd.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}



// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
int Proc1FA(string& outname, string& inpath, int imin, int imax)
{
  vector<string> infiles;
  char nbuff[1024];
  for(int ii=imin; ii<=imax; ii++) {
    sprintf(nbuff,"%s/signal%d.fits",inpath.c_str(),ii);
    infiles.push_back(nbuff);
  }
  uint_4 nmaxz;
  uint_4 paqsz, npaqf;
  if (fgrdfits) {
    DecodeMiniFitsHeader(infiles[0],paqsz, npaqf, fgnotrl); 
    nmaxz = infiles.size()*npaqf/NPaqinZone;
    cout << " mcrd/Proc1FRawA/ReadFits: NZones=" << NZones << " NbPaq=" << NPaqinZone 
         << " NPaq in file=" << npaqf
         << " PaqSz=" << paqsz << " NMaxZ=" << nmaxz << " NMean=" << NMean << endl;
    if (fgraw) cout << " ======>  RAW data processing " << endl;
    else cout << " ======>  FFT data processing " << endl;
  }
  else {
    paqsz = GPaqSz;
    nmaxz = NGenZ;
    cout << " mcrd/Proc1FRawA/GeneratePaquets: NZones=" << NZones << " NbPaq=" << NPaqinZone 
         << "  PaqSz=" << paqsz << " NMean=" << NMean << endl;
  }
  RAcqMemZoneMgr mmgr(NZones, 1, NPaqinZone, paqsz, sizeof(complex<r_4>)*paqsz);
  mmgr.SetFinalizedMask((uint_4)MemZS_ProcA);

  BRFitsReader reader(mmgr, infiles, fgnotrl);
  TestPCIWrapperNODMA pciw(paqsz,LossRate);
  PCIEReader pcird(pciw, paqsz, paqsz, mmgr, nmaxz, BR_Copy);

  outname += "/Ch12";
  BRProcA2C proc(mmgr, outname, fgraw, NMean, nmaxz, fghist, NFreqSMap, fgnotrl);

  cout << " mcrd/Proc1FA: Starting threads (reader, proc) ... " << endl;

  if (fgrdfits) reader.start();
  else pcird.start();

  proc.start();
  sleep(1);
  cout << " mcrd/Proc1FA: Waiting for reader thread to finish ... " << endl;
  if (fgrdfits) reader.join();
  else pcird.join();
  cout << " mcrd/Proc1FA: Reader finished, waiting for process thread to finish ... " << endl;
  sleep(2);
  mmgr.Stop();
  proc.join();
  mmgr.Print(cout);
  return 0;
}


// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
int Proc2FAB(string& outname, string& path1, string& path2, int imin, int imax) 
{
  vector<string> infiles1;
  vector<string> infiles2;
  char nbuff[1024];
  for(int ii=imin; ii<=imax; ii++) {
    sprintf(nbuff,"%s/signal%d.fits",path1.c_str(),ii);
    infiles1.push_back(nbuff);
    sprintf(nbuff,"%s/signal%d.fits",path2.c_str(),ii);
    infiles2.push_back(nbuff);
  }
  uint_4 nmaxz;
  uint_4 paqsz, npaqf;
  if (fgrdfits) {
    DecodeMiniFitsHeader(infiles1[0],paqsz, npaqf, fgnotrl); 
    nmaxz = infiles1.size()*npaqf/NPaqinZone;
    cout << " mcrd/Proc2FAB/ReadFits: NZones=" << NZones << " NbPaq=" << NPaqinZone 
         << " NPaq in file=" << npaqf
         << " PaqSz=" << paqsz << " NMaxZ=" << nmaxz << " NMean=" << NMean << endl;
    if (fgraw) cout << " ======>  RAW data processing " << endl;
    else cout << " ======>  FFT data processing " << endl;
  }
  else {
    paqsz = GPaqSz;
    nmaxz = NGenZ;
    cout << " mcrd/Proc2FAB/GeneratePaquets: NZones=" << NZones << " NbPaq=" << NPaqinZone 
         << "  PaqSz=" << paqsz << " NMean=" << NMean << endl;
  }
  RAcqMemZoneMgr mmgr1(NZones, 1, NPaqinZone, paqsz, sizeof(complex<r_4>)*paqsz);
  mmgr1.SetFinalizedMask((uint_4)(MemZS_ProcA)|(uint_4)(MemZS_ProcB));
//  mmgr1.SetFinalizedMask((uint_4)(MemZS_ProcA));
  RAcqMemZoneMgr mmgr2(NZones, 1, NPaqinZone, paqsz, sizeof(complex<r_4>)*paqsz);
  mmgr2.SetFinalizedMask((uint_4)(MemZS_ProcA)|(uint_4)(MemZS_ProcB));
//  mmgr2.SetFinalizedMask((uint_4)(MemZS_ProcA));

//-- Lecture des fichiers
  BRFitsReader reader1(mmgr1, infiles1, fgnotrl);
  BRFitsReader reader2(mmgr2, infiles2, fgnotrl);

//-- Generation de paquets 
  TestPCIWrapperNODMA pciw1(paqsz,LossRate);
  PCIEReader pcird1(pciw1, paqsz, paqsz, mmgr1, nmaxz, BR_Copy);
  TestPCIWrapperNODMA pciw2(paqsz,LossRate);
  PCIEReader pcird2(pciw2, paqsz, paqsz, mmgr2, nmaxz, BR_Copy);

  string outname1 = outname;
  outname1 += "/Ch12";
  BRProcA2C proc1(mmgr1, outname1, fgraw, NMean, nmaxz, fghist, NFreqSMap, fgnotrl);
  string outname2 = outname;
  outname2 += "/Ch34";
  BRProcA2C proc2(mmgr2, outname2, fgraw, NMean, nmaxz, fghist, NFreqSMap, fgnotrl,2);
  string outname12 = outname;
  outname12 += "/Ch1234";
  BRProcB4C proc12(mmgr1, mmgr2, outname12, fgraw, NMean, nmaxz, fgnotrl);

  cout << " mcrd/Proc2FRawAB: Starting threads (reader1,2, procA1,2, procAB) ... " << endl;
//  cout << "[1]--- CR to continue ..." << endl;   char ans[32]; gets(ans); 
  if (fgrdfits) { 
    reader1.start();
    reader2.start();
  }
  else { 
    pcird1.start();
    pcird2.start();
  }
//  sleep(1);  mmgr1.Print(cout);  mmgr2.Print(cout);
  proc1.start();
  proc2.start();
  proc12.start();

  sleep(1);
  cout << " mcrd/Proc2FRawAB: Waiting for reader threads to finish ... " << endl;
  if (fgrdfits) { 
    reader1.join();
    reader2.join();
  }
  else { 
    pcird1.join();
    pcird2.join();
  }
  cout << " mcrd/Proc2FRawAB: Readers finished, waiting for process thread to finish ... " << endl;
  sleep(2);
  mmgr1.Stop();
  mmgr2.Stop();
  proc1.join();
  proc2.join();
  proc12.join();
  mmgr1.Print(cout);
  mmgr2.Print(cout);
  return 0;
}


/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- mcrd.cc : Reading/Processing BAORadio FITS files" << endl;
  cout << " Usage:  mcrd ACT OutPath InPath [InPath2] Imin,Imax\n" 
       << "              [NZones,NPaqinZone] [NMean] [NFreqSMap]\n"
       << "              [GPaqSz,NGenZones,LossRate]" << endl;
  if (fgshort) {
    cout << " mcrd -h for detailed instructions" << endl;
    return 1;
  }
  cout << " ACT= -2[ghnf] -> 1 fiber, 2 channel processing (ProcA)\n"
       << " ACT= -4[ghnf] -> 2 fibers, 4 channels (ProcA, ProcB)\n"
       << "   f FFT data ( raw if not) -> \n" 
       << "   n (notrl) -> FITS files without frame trailer \n" 
       << "   g (generate paquets) -> generate paquets instead of reading fits files\n"
       << "   h (time sample histograms) -> compute time sample histograms also \n" 
       << "   Example: ACT = -2h  1 fiber, 2 raw channels and compute time sample histograms" <<endl; 
  cout << " OutPath : Output directory name " << endl;
  cout << " InPath [InPath2] Imin,Imax: Input fits files directory name(s)\n"
       << "    (Inpath=Fiber1, InPath2=Fiber2) and sequence numbers \n" 
       << "    FileNames=InPath/signalII.fits Imin<=II<=Imax \n"
       << " NZones,NbPaqinZone : Number of Zones and number of paquets in one zone\n"
       << "     of the RAcqMemZoneMgr memory manager (default = 4,128)\n"
       << " NMean: Number of packet used for spectra/visibility computation (def=1024)\n"
       << " NFreqSMap: Frequency binning for 2D (time-freq) spectra maps (def=0->No maps)\n"
       << " GPaqSz,NGenZones,LossRate: Paquet Size, Number of memory zones filled and\n"
       << "      loss rate (-2g,-4g) , default=16424,100,0.1" << endl;
  return 1;
}