source: Sophya/trunk/AddOn/TAcq/visfits2dt.cc@ 4060

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

/* ---------------------------------------------------------- 
   Projet BAORadio - (C) LAL/IRFU  2011

   Programme de lecture des fichiers matrices de visibilites 
   au format fits pour ecrire une DataTable de visib_moyenne(temps)
   Fichier de visibilites produits par le prog vismf 
   R. Ansari, C. Magneville   -  LAL/Irfu
   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 <utilarr.h>
#include "fitsarrhand.h"
#include "fitshdtable.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 ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
                vector<sa_size_t> tfrlist);
int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
                  sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;

//------------------------------------------------------------------------------------------------------------

/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- visfits2dt : Read fits visibility matrices produced by visfmib to make \n"
       << "     time-frequency matrices and Visibilites=f(time) datatable " << endl;
  cout << " Usage: visfits2dt InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList [DT] \n" << endl;
  if (fgshort) {
    cout << " svv2mtx -h for detailed instructions" << endl;
    return 1;
  }
  cout << " InOutPath : Input/Output directory name \n" 
       << " Imin,Imax,IStep: Input FITS  files sequence number (vismtxII.fits) \n"
       << "    FileNames=InOutPath/vismtxII.fits Imin<=II<=Imax II+=IStep \n" 
       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for DataTable\n"
       << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n" 
       << " DT : Fill datatable if DT arg specified " << endl;
  return 1;
}

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

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

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

}

//-------------------------------------------------------------------- 
//       Traitement fichiers produits par vismfib  (V Nov09)
/* --Fonction-- */
static  bool fgfilldt=false;

int DecodeProc(int narg, char* arg[])
{
  // Decodage des arguments et traitement 
  string inoutpath = arg[0];    
  int imin=0;
  int imax=0;
  int istep=1;
  sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
  int jf1=0;
  int jf2=0;
  int nfreq=0;
  sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
  int card=1;
  vector<sa_size_t> rowlist;
  EnumeratedSequence es;
  int nbad;
  es.Append(arg[3], nbad, ",");
  for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
  fgfilldt=false;
  if ((narg>4)&&(strcmp(arg[4],"DT")==0))  fgfilldt=true;
  //  if (rowlist.size()<1) rowlist.push_back(0);

  cout << " ----- visfits2dt/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
       << imin << "," << imax << "," << istep << " Card=" << card << endl;
  cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
  cout << " RowList: " ;
  for(int j=0; j<rowlist.size(); j++)  cout << rowlist[j] << " , " ;
  cout << endl;
  int rc=ProcVisMtxFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
  return rc;
}

// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
/* --Fonction-- */
int ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
                vector<sa_size_t> rowlist)
{
  Timer tm("ProcVisMtxFiles");
  
  DataTable visdt;
  visdt.AddDoubleColumn("mfc");
  visdt.AddDoubleColumn("mtt");
  visdt.AddIntegerColumn("jfreq");
  visdt.AddIntegerColumn("numch");
  visdt.AddFloatColumn("vre");
  visdt.AddFloatColumn("vim");

  vector< TMatrix< r_4 > > vmtf;
  vector< TMatrix< r_4 > > vmtfa;

  if (jf1<1)  jf1=1;
  if ((jf2<1)||(jf2<jf1))  jf2=jf1; 
  if (nfreq<1) nfreq=1;
  int djf=(jf2-jf1)/nfreq;
  if (djf<1) djf=1;

  cout << " --- ProcVisMtxFiles Frequency binng Start=" << jf1 << " End= " << jf2 << " DFreq=" << djf << " NBinFreq=" << nfreq << endl;
  char fname[1024];

  cout << " ---- ProcVisMtxFiles()-Begin - Reading chanum vector" << endl;
  TVector< int_4 > chanum;      
  {
  sprintf(fname, "%s/chanum.fits",inoutpath.c_str());
  FitsInOutFile fic(fname, FitsInOutFile::Fits_RO);
  fic.MoveAbsToHDU(3);
  fic >> chanum;
  }
  
  cout << " ---- ProcVisMtxFiles()-Read chanum done " << endl;
  sa_size_t nrows = (imax-imin+1)/istep;
  sa_size_t ktime=0;

  //  sa_size_t mtf_binfreq=25;
  sa_size_t mtf_bintime=1;

  sa_size_t ncols=1;

  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    sprintf(fname, "%s/vismtx%d.fits",inoutpath.c_str(),ifile);
    cout << " ProcVisMtxFiles[" << ifile << "] opening file " << fname << endl;
    FitsInOutFile fin(fname, FitsInOutFile::Fits_RO);
    TMatrix< r_4 > vismtx;

    fin >> vismtx;
    //    cout << vismtx.Info();

    if ((ifile==imin)&&(rowlist.size()>0)) {
      ncols = vismtx.NCols();
      cout << " ProcVisMtxFiles/Info: Output Time-Frequency matrices NRows(time) " 
           << nrows/mtf_bintime+1 << " NCols=2*NFreq) =" << nfreq << endl;
      for(size_t j=0; j<rowlist.size(); j++) {
        vmtf.push_back(TMatrix< r_4 >(nrows/mtf_bintime+1, nfreq*2));
        vmtfa.push_back(TMatrix< r_4 >(nrows/mtf_bintime+1, nfreq));
      }
    }

    if (rowlist.size()>0) {
      for(size_t j=0; j<rowlist.size(); j++) {
        sa_size_t jfreb=0;
        for(sa_size_t jf=jf1; jf<jf2; jf++) {
          sa_size_t jfreb=(jf-jf1)/djf;
          if (jfreb>=nfreq)  break;
          vmtf[j](ktime/mtf_bintime,jfreb*2)+=vismtx(rowlist[j],jf*2); // partie reelle 
          vmtf[j](ktime/mtf_bintime,jfreb*2+1)+=vismtx(rowlist[j],jf*2+1); // partie imaginaire 
          vmtfa[j](ktime/mtf_bintime,jfreb)+=vismtx(rowlist[j],jf*2)*vismtx(rowlist[j],jf*2); // module carre - z.real^2 
          vmtfa[j](ktime/mtf_bintime,jfreb)+=vismtx(rowlist[j],jf*2+1)*vismtx(rowlist[j],jf*2+1); // module carre z.imag^2 
        }
      }
    }
    //    cout << " DBG* kr=" << kr << " kr/mtf_bintime=" << kr/mtf_bintime 
    //     << " ncols/2/mtf_binfreq=" << ncols/2/mtf_binfreq << endl;
    ktime++;

     
// Calcul moyenne dans des bandes en frequence 
    if (fgfilldt)   FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
  }

  if (rowlist.size()>0) {
    sprintf(fname, "!%s/vistfreqmtx.fits",inoutpath.c_str());
    cout << "ProcVisMtxFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;  
    FitsInOutFile fo(fname, FitsInOutFile::Fits_Create);
    for(size_t j=0; j<rowlist.size(); j++)  {
      cout << " Writing Time-Freq visibility matrix[" << j << "] for VisIJ= " << chanum(rowlist[j]) << endl;
      fo << vmtf[j];
    }
    sprintf(fname, "!%s/acortfreqmtx.fits",inoutpath.c_str());
    cout << "ProcVisMtxFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;  
    FitsInOutFile foa(fname, FitsInOutFile::Fits_Create);
    for(size_t j=0; j<rowlist.size(); j++)  {
      cout << " Writing Time-Freq modulesquare matrix[" << j << "] for Mod2IJ= " << chanum(rowlist[j]) << endl;
      foa << vmtfa[j];
    }

  }
  if (fgfilldt)  {
    cout << visdt;
    sprintf(fname, "!%s/visidt.fits",inoutpath.c_str());
    cout << "ProcVisMtxFiles: writing visibility datatable to file " << fname << endl;
    FitsInOutFile fod(fname, FitsInOutFile::Fits_Create);
    fod << visdt;
  }
  return 0;
}

/* --Fonction-- */
int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
                  sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) 
{
  double xnt_[20];
  xnt_[0]=(double)vismtx_.Info()["MEANFC"];  
  xnt_[1]=(double)vismtx_.Info()["MEANTT"]/1.25e8;

  uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
  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*2)/(r_4)(nmean_);
        xnt_[5]=vismtx_(rv,jf*2+1)/(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*2);
          moyimag+=vismtx_(rv,jjf*2+1);
        }
        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;
}