Added msvis2dt.cc to read meanspecII.ppf files (mfacq with monitroring...

Added msvis2dt.cc to read meanspecII.ppf files (mfacq with monitroring activated) and vismtxII.ppf files - DEBUG configuration of acquisition with 1 ADC board with RAW firware and one ADC board with FFT firmware, Reza 25/9/2015

msvis2dt.cc

+// Utilisation de SOPHYA pour faciliter les tests ...
+#include "sopnamsp.h"
+#include "machdefs.h"
+
+/* ---------------------------------------------------------- 
+   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2015 
+
+   msvis2dt: programme de lecture des fichiers meanspecII.ppf 
+   (sortie de monitoring), et vismtxII.ppf (visibilities) 
+   Configuration de debug avec une carte avec firmware RAW et 
+   une avec firmware FFT - Septembre 2015 
+   R. Ansari, J.E. Campagne, C. Magneville   -  LAL/Irfu
+   ---------------------------------------------------------- */
+
+// 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 sophya mesure ressource CPU/memoire ...
+#include "resusage.h"
+#include "ctimer.h"
+#include "timing.h"
+
+// include lecteur de fichiers visibilites 
+// #include "visip4reader.h"
+
+//--------------------------- Fonctions de ce fichier   ------------------- 
+int Usage(bool fgshort=true);
+// int DecodeArgs(int narg, char* arg[]);
+
+/* --Fonction-- */
+int Usage(bool fgea)
+{
+  cout << " --- msvis2dt.cc : Read PPF files produced by mfacq (meanspecII.ppf vismtxII.ppf) \n" << endl;
+  if (fgea) cout << " msvis2dt: Missing or wrong argument ! " << endl;
+  cout << " Usage: msvis2dt [-tfm Nf] InPath_vismtx InPath_meanspec Imin,Imax OutPPFFile Freq1 Freq2 DeltaFreq [PrtLev=0] \n"
+       << "  -tfm Nf : create time-frequncy maps of auto-correlation binning Nf frequencies \n" 
+       << "  InPath_vismtx: path for input vismtxII.ppf (on bao5) \n"
+       << "  InPath_meanspec: path for input meanspecII.ppf (on bao6) \n"
+       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
+       << "  OutPPFFile: Output PPF file name \n"
+       << "  Freq1 Freq2: the two frequencies for which cross-correlations are computed [=1375 1410 MHz] \n"
+       << "  DeltaFreq: frequency band (in MHz) [=1 MHz] \n"
+       << "    Note: bande_freq=[freq1,2-deltafreq , freq1,2+deltafreq] \n"<< endl;
+  /*
+  if (fgshort) {
+    cout << " rdvisip4 -h for detailed instructions" << endl;
+    return 1;
+  }
+  */
+  return 1;
+}
+
+//----------------------------------------------------
+//----------------------------------------------------
+int main(int narg, char* arg[])
+{
+  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
+  if (narg<9) return Usage(true);
+  int offa=0;
+  bool FgTFM=false;   // true -> create time-frequency maps of cross-correlations 
+  sa_size_t TFMfbin=8;
+  if (strcmp(arg[1],"-tfm")==0) {
+    if (narg<11) return Usage(true);
+    FgTFM=true;
+    TFMfbin=atoi(arg[2]);
+    if (TFMfbin<1) TFMfbin=8;
+    offa=2;
+  }
+  string path5 = arg[offa+1];
+  string path6 = arg[offa+2];
+  int Imin,Imax,Istep=1; 
+  sscanf(arg[offa+3],"%d,%d",&Imin,&Imax);
+  Istep=1;
+  string outfile=arg[offa+4];
+  // frequency range definition 
+  double deltafreq=1.;
+  double freq1=1375.;
+  double freq2=1410.;
+  int nargo=narg-offa;
+  if ((nargo>5)&&(strcmp(arg[offa+5],"-")!=0))  freq1=atof(arg[offa+5]);
+  if ((nargo>6)&&(strcmp(arg[offa+6],"-")!=0))  freq2=atof(arg[offa+6]);
+  if ((nargo>7)&&(strcmp(arg[offa+7],"-")!=0))  deltafreq=atof(arg[offa+7]);
+  // time averaging interval definition, in term of visibility files
+  int deltaIavg=1;
+  //  if ((nargo>8)&&(strcmp(arg[offa+8],"-")!=0))  deltaIavg=atoi(arg[offa+8]);
+  if (deltaIavg<1)  deltaIavg=1;
+  int prtlev=1;
+  if (nargo>8) prtlev=atoi(arg[offa+8]);
+
+  cout << " msvis2dt/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
+       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
+       << " freq1="<<freq1<<" freq2="<<freq2<<" deltafreq="<<deltafreq<<" MHz\n"
+    //       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
+       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
+
+  string TFMoutfile;
+  if (FgTFM) {
+    size_t l = outfile.length();
+    size_t pp = outfile.rfind('.');
+    if (pp<l) TFMoutfile = outfile.substr(0,pp) + "_tfm.ppf";
+    cout << " msvis2dt/Info: Time-Frequency maps of Xcor will be created and saved to file "
+	 <<TFMoutfile<<endl;
+  }
+
+  //---- calcul des index de bandes de frequences 
+  // pas en frequence du firmware FFT 
+  double deltanufft=250./4096;    // 250 MHz en 4096 frequences 
+  sa_size_t DJF=deltafreq/deltanufft;
+  double freqstart=1250.;   // Bande de 1250-1500 MHz 
+  sa_size_t JFmin=(freq1-1250.-deltafreq)/deltanufft;
+  sa_size_t JFmax=(freq1-1250.+deltafreq)/deltanufft;
+  sa_size_t JFminB=(freq2-1250.-deltafreq)/deltanufft;
+  sa_size_t JFmaxB=(freq2-1250.+deltafreq)/deltanufft;
+  cout << " msvis2dt/Info:   frequency bands ... " <<endl;
+  cout << " Band 1: "<<freq1<<" +/-"<<deltafreq<<" MHz -> "<<JFmin<<" <=NumFreq<= "<<JFmax<<endl;
+  cout << " Band 2: "<<freq2<<" +/-"<<deltafreq<<" MHz -> "<<JFminB<<" <=NumFreq<= "<<JFmaxB<<endl;
+
+  // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792) 
+  sa_size_t JFmin21=2776,JFmax21=2808;   // +/- 1 MHz autour de 1420 MHz 
+  sa_size_t JFmin21_5=2710,JFmax21_5=2874;  //  +/- 5 MHz autour de 1420 MHz 
+  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
+       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
+ 
+  HiStatsInitiator _inia;
+  //   TArrayInitiator  _inia;
+
+  int rc = 0;
+  try {
+    ResourceUsage resu;
+
+    //  Numero des lignes des auto-correlations dans la matrice des visibilites 
+    sa_size_t KVAC[4]={0,4,7,9}; 
+    //  Numero des lignes des cross-correlations 1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H dans la matrice des visibilites 
+    sa_size_t KVCXHH[6]={1,2,3,5,6,8}; 
+
+    const char* dtnames[43]={"k","day","time",
+			     "p1h","p2h","p3h","p4h",
+			     "p1hB","p2hB","p3hB","p4hB",
+			     "p1hHI","p2hHI","p3hHI","p4hHI",
+			     "p1hHI5","p2hHI5","p3hHI5","p4hHI5",
+                             "xh1h2","xh1h3","xh1h4","xh2h3","xh2h4","xh3h4",
+                             "xh1h2B","xh1h3B","xh1h4B","xh2h3B","xh2h4B","xh3h4B",
+			     "xh1h2HI","xh1h3HI","xh1h4HI","xh2h3HI","xh2h4HI","xh3h4HI",
+			     "dayraw","timeraw",
+			     "p3hraw","p3hrawB","p3hrawHI","p3hrawHI5"};
+
+    DataTable dt(64);
+    dt.AddIntegerColumn(dtnames[0]);
+    dt.AddIntegerColumn(dtnames[1]);
+    dt.AddFloatColumn(dtnames[2]);
+    for(int kk=0; kk<16; kk++) 
+      dt.AddFloatColumn(dtnames[3+kk]);
+    for(int kk=0; kk<18; kk++) 
+      dt.AddComplexColumn(dtnames[3+16+kk]);
+
+    dt.AddIntegerColumn(dtnames[3+16+18]);    
+    dt.AddFloatColumn(dtnames[3+16+18+1]);
+    for(int kk=0; kk<4; kk++) 
+      dt.AddFloatColumn(dtnames[3+16+18+2+kk]);
+
+    Range freqs(JFmin,JFmax);
+    Range freqs21(JFmin21,JFmax21);
+    Range freqs21_5(JFmin21_5,JFmax21_5);
+    Range freqsB(JFminB,JFmaxB);
+
+    
+    double acrdt[4];  // 4 bandes x 1 autocor raw/meanspec 
+    double acdt[16];   // 4 bandes x [ 4 autocor vismtx ]
+    complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 
+
+    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
+    double mttag;
+    int cnt=0;
+    //----- 5 TimeFrequency maps of auto-correlations 
+    vector< TArray< r_4 > > vtfm;
+    
+    // Getting en empty row_ptr
+    DataTableRowPtr rowp = dt.EmptyRowPtr();
+
+    //---- for the time-freqency map filling    
+    sa_size_t TFMtmidx=0;
+    sa_size_t tfmSX, tfmSY;
+
+    char filename[1024];
+    for(int jf=Imin; jf<=Imax; jf+=Istep) {
+      sprintf(filename,"%s/vismtx%d.ppf",path5.c_str(),jf);
+      if (prtlev>1)  cout << "JFile="<<jf<<" -> Opening file " << filename << endl;
+      PInPersist pin5(filename);
+      TMatrix< complex<float> >  vismtx;
+      time_t tm5 = pin5.CreationDate();
+      struct tm * ptms5 = gmtime(&tm5);
+      TimeStamp tms5(ptms5->tm_year, ptms5->tm_mon, ptms5->tm_mday, ptms5->tm_hour, ptms5->tm_min, ptms5->tm_sec);
+      pin5 >> vismtx;
+
+      sprintf(filename,"%s/meanspec%d.ppf",path6.c_str(),jf);
+      if (prtlev>1)  cout << "JFile="<<jf<<" -> Opening file " << filename << endl;
+      PInPersist pin6(filename);
+      TMatrix<float> spectre;
+      time_t tm6 = pin6.CreationDate();
+      struct tm * ptms6 = gmtime(&tm6);
+      TimeStamp tms6(ptms6->tm_year, ptms6->tm_mon, ptms6->tm_mday, ptms6->tm_hour, ptms6->tm_min, ptms6->tm_sec);
+      pin6 >> PPFNameTag("spectre") >> spectre;
+
+      if (vismtx.NCols() != spectre.NCols())  
+	throw PError("msvis2dt/ERROR: vismtx and spectre do NOT have same number of frequencies");
+      if (jf == Imin) { 
+	if (FgTFM) {  //allocating time-frequency maps 
+	  tfmSX=(Imax-Imin+1)/Istep/deltaIavg;
+	  tfmSY=vismtx.NCols()/TFMfbin;
+	  cout << " msvis2dt/Info: allocating Time-Freqncy maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
+	  for(int k=0; k<5; k++) vtfm.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
+	}
+      }
+      //---------------------------------------------------------------------------------
+      // pour les 4 autocorrelations provenant des visibilites  
+      for(int k=0; k<4; k++)  {
+	TVector< complex<float> > acv = vismtx.Row( KVAC[k] );
+	acdt[k] = acv(freqs).Sum().real();       // integration en bande 1
+	acdt[k+4] = acv(freqsB).Sum().real();      // integration en bande large B
+	acdt[k+8] = acv(freqs21).Sum().real();     // integration a +/- 1 MHz @ 1420 MHz 
+	acdt[k+12] = acv(freqs21_5).Sum().real();   // integration a +/- 5 MHz @ 1420 MHz 
+
+	if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
+	  TArray< r_4 > & tfmap = vtfm[k];
+	  for(sa_size_t jy=0; jy<tfmSY; jy++) {
+	    complex<r_4> zz = acv( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
+	    tfmap(TFMtmidx, jy) = zz.real();
+	  } 
+	}  //----- end of time-frequency maps filling 
+      }
+      {
+      // Voie 3H-raw 
+	TVector<float> acr = spectre.Row(0);
+	acrdt[0] = acr(freqs).Sum();    // integration en bande 1
+	acrdt[1] = acr(freqsB).Sum();    // integration en bande large B
+	acrdt[2] = acr(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
+	acrdt[3] = acr(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
+	if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
+	  TArray< r_4 > & tfmap = vtfm[4];
+	  for(sa_size_t jy=0; jy<tfmSY; jy++) {
+	    complex<r_4> zz = acr( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
+	    tfmap(TFMtmidx, jy) = zz.real();
+	  } 
+	}  //----- end of time-frequency maps filling 
+      }
+      //---- pour les 6 cross-cor
+      for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
+	TVector< complex<float> > cxv = vismtx.Row( KVCXHH[k] );
+	cxdt[k] = cxv(freqs).Sum();    // integration en bande 1
+	cxdt[k+6] = cxv(freqsB).Sum();    // integration en bande large B
+	cxdt[k+12] = cxv(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
+      }
+      //  On calcule les moyennes par bande de frequence
+      // Moyennes dans la bande 1
+      //       double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
+      double fnorm=1./(double)(JFmax-JFmin+1);
+      for(int k=0; k<4; k++)  acdt[k]*=fnorm; 
+      acrdt[0]*=fnorm; 
+      for(int k=0; k<6; k++)  cxdt[k]*=complex<r_8>(fnorm,0.);  
+      // Moyennes dans la bande 2 
+      fnorm=1./(double)(JFmaxB-JFminB+1);
+      for(int k=0; k<4; k++)  acdt[k+4]*=fnorm; 
+      acrdt[1]*=fnorm; 
+      for(int k=0; k<6; k++)  cxdt[k+6]*=complex<r_8>(fnorm,0.);   
+      // Moyennes ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420 
+      fnorm=1./(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
+      for(int k=0; k<8; k++)  acdt[k+12] = (acdt[k+12]-acdt[k+8])*fnorm;
+      acrdt[3] = (acrdt[3]-acrdt[2])*fnorm; 
+      // Moyennes ds +/1 MHz @1420 
+      fnorm=1./(double)(JFmax21-JFmin21+1);
+      for(int k=0; k<4; k++)  acdt[k+8]*=fnorm; 
+      acrdt[2]*=fnorm; 
+      for(int k=0; k<6; k++)  cxdt[k+12]*=complex<r_8>(fnorm,0.);   
+      
+      // remplissage datatable 
+      dt.NextRowPtr(rowp);
+      rowp(0)=(int_4)cnt;   
+      rowp(1)=(int_4)(tms5.DaysPart()-dateorg.DaysPart());   
+      rowp(2)=(r_4)tms5.SecondsPart();
+      for(int k=0; k<16; k++)   rowp(3+k) = (float)acdt[k];
+      for(int k=0; k<18; k++)   rowp(19+k) = complex<float>((float)cxdt[k].real(), (float)cxdt[k].imag());
+      rowp(37)=(int_4)(tms6.DaysPart()-dateorg.DaysPart());   
+      rowp(38)=(r_4)tms6.SecondsPart();
+      for(int k=0; k<4; k++)   rowp(39+k) = (float)acrdt[k];
+      //  ... done 
+      cnt++;
+      
+      if (prtlev>1) cout << "JFile="<<jf<<" DateFile5:"<<tms5<<" DateFile6:"<<tms6<<" p3hraw="<<acrdt[0]<<endl;
+    }
+ 
+    cout << " msvis2dt/Info: count="<<cnt<<" visimtx/meanspec read "<<endl;
+    dt.SetShowMinMaxFlag(true);
+    cout << dt;
+    cout<<"  msvis2dt/Info: Saving auto/cross-corr=f(time) datatable to PPF file "<<outfile<<endl;
+    POutPersist po(outfile);
+    po<<dt; 
+
+    if (FgTFM) {   // --- renormalizing and saving time-frequency maps 
+      cout<<"  msvis2dt/Info: Saving time-frequency maps to PPF file "<<TFMoutfile<<endl;
+      POutPersist potfm(TFMoutfile);
+      const char* tfm_names[6]={"TFM_3Hraw", "TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H"};
+      for(int k=0; k<5; k++)  {   // loop over the 6 Xcor 
+	TArray< r_4 > & tfmap = vtfm[k];
+	tfmap /= (r_4)(1./(double)deltaIavg/(double)TFMfbin);
+	potfm << PPFNameTag(tfm_names[k]) << tfmap;
+      }
+    }
+    //    resu.Update();
+    cout << resu;   // Update est fait lors du print
+  }
+  catch (PException& exc) {
+    cerr << " msvis2dt.cc catched PException " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n msvisi2dt.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+    rc = 78;
+  }
+  catch (...) {
+    cerr << " msvisi2dt.cc catched unknown (...) exception  " << endl; 
+    rc = 79; 
+  } 
+
+  cout << ">>>> msvisi2dt.cc ------- END ----------- RC=" << rc << endl;
+  return rc;
+
+}
+
+