Changeset 3701 in Sophya for trunk/AddOn

Timestamp:

Dec 5, 2009, 1:00:11 AM (16 years ago)

Author:

ansari

Message:

MAJ du prog svv2mtx pour traitement fichiers produits par vismfib, Reza 4/12/09

Location:

trunk/AddOn/TAcq

Files:

• brproc.cc (modified) (4 diffs)
• svv2mtx.cc (modified) (8 diffs)

trunk/AddOn/TAcq/brproc.cc

@@ -57 +57 @@
   visdt_.AddFloatColumn("vim");
 
-  /*
-  char colnm[32];
-  for(size_t i=0; i<2*memgr_.NbFibres(); i++) {
-    for(size_t j=i; j<2*memgr_.NbFibres(); j++) {
-      sprintf(colnm,"vr_%d_%d",(int)(i+1),(int)j+1);
-      visdt_.AddFloatColumn(colnm);
-      sprintf(colnm,"vi_%d_%d",(int)(i+1),(int)j+1);
-      visdt_.AddFloatColumn(colnm);
-    }
-  }
-  */
   if (nmean_ < 1)  nmean_=memgr_.NbPaquets();
   if (nmean_ < 1)  nmean_=1;
@@ -95 +84 @@
   POutPersist po(outpath_+"visdt.ppf");
   po << visdt_;
+  POutPersist poc(outpath_+"chanum.ppf");
+  po << chanum_;
+
   delete[] xnt_;
 }
@@ -109 +101 @@
   if (totnbpaq_%nmean_ == 0) {
     if (totnbpaq_ > 0) {
+      moyfc_/=nmean_;
+      moytt_/=nmean_;
+      vismtx_.Info()["MeanFC"] = moyfc_;
+      vismtx_.Info()["MeanTT"] = moytt_;
+      vismtx_.Info()["NPAQSUM"] = nmean_;
+
       // ATTENTION : Matrice visibilites non moyennee 
       char nfile[32];
@@ -117 +115 @@
       cout << numfile_ << "-BRVisibilityCalculator::Process() NPaqProc="
            << totnbpaq_ << "  -> Visibility Matrix in " << flnm << endl;
-      moyfc_/=nmean_;
-      moytt_/=nmean_;
       FillVisibTable(moyfc_, moytt_);
       numfile_++;

trunk/AddOn/TAcq/svv2mtx.cc

@@ -4 +4 @@
 
 /* ---------------------------------------------------------- 
-   Programme de lecture multi canaux pour BAORadio  
-   R. Ansari, C. Magneville
+   Projet BAORadio - (C) LAL/IRFU  2008-2010
+
+   Programme de lecture des fichiers vecteurs/matrices de 
+   visibilites produits par mcrd / vismfib 
+   R. Ansari, C. Magneville   -  LAL/Irfu
    V : Mai 2009
    ---------------------------------------------------------- */
@@ -23 +26 @@
 // #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 ...
@@ -33 +38 @@
 
 
+//--------------------------- Fonctions de ce fichier   ------------------- 
 int Usage(bool fgshort=true);
-// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
-int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card=1);
-int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq);
+
+int DecodeProc(int narg, char* arg[]);
+int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
+                vector<sa_size_t> tfrlist);
+int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > visdt_, TVector< uint_4> chanum_,
+                  sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
+
+int DecodeProcVJun09(int narg, char* arg[]);
+int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card);
+//------------------------------------------------------------------------------------------------------------
+
+/* --Fonction-- */
+int Usage(bool fgshort)
+{
+  cout << " --- svv2mtx.cc : Read PPF files produced by mcrd/visfmib to make time-frequency\n"
+       << "     matrices and Visibilites=f(time) datatable " << endl;
+  cout << " Usage: svv2mtx -mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1] \n" 
+       << "   OR \n" 
+       << "        svv2mtx InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList" << endl; 
+  if (fgshort) {
+    cout << " svv2mtx -h for detailed instructions" << endl;
+    return 1;
+  }
+  cout << " -mcrd : Read mcrd output files -> ProcSVFilesVJun09() \n" 
+       << " InOutPath : Input/Output directory name \n" 
+       << " Imin,Imax,IStep: Input PPF files sequence number \n"
+       << "    FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n" 
+       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
+       << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n"
+       << " card=1 Ch12 , card=2 Ch34    (mcrd output files) " << endl;
+  return 1;
+}
 
 //----------------------------------------------------
@@ -44 +79 @@
   if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
   if (narg<4) return Usage(true);
-
   HiStatsInitiator _inia;
   //   TArrayInitiator  _inia;
@@ -50 +84 @@
   int rc = 0;
   try {
-    string inoutpath = 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;
-    sscanf(arg[3],"%d,%d,%d",&jf1,&jf2,&nfreq);
-    int card=1;
-    if (narg>4) card=atoi(arg[4]); 
-    cout << " ----- svv2mtx.cc Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
-         << imin << "," << imax << "," << istep << " Card=" << card << endl;
-    cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
+    bool fgvjun09=false;
+    if ((strcmp(arg[1],"-mcrd")==0))  fgvjun09=true;
     ResourceUsage resu;
-    rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
+    if (fgvjun09)  DecodeProcVJun09(narg-2, arg+2);
+    else  DecodeProc(narg-1, arg+1);
     resu.Update();
     cout << resu;
@@ -89 +112 @@
 }
 
-
+//-------------------------------------------------------------------- 
+//       Traitement fichiers produits par vismfib  (V Nov09)
+/* --Fonction-- */
+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;
+  if (narg>3) { 
+    EnumeratedSequence es;
+    int nbad;
+    es.Append(arg[3], nbad, ",");
+    for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
+  }
+  if (rowlist.size()<1) rowlist.push_back(0);
+
+  cout << " ----- svv2mtx/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=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
+  return rc;
+}
 
 // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
+/* --Fonction-- */
+int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
+                vector<sa_size_t> rowlist)
+{
+  Timer tm("ProcSVFiles");
+
+  vector< TMatrix< complex<r_4> > > vmtf(rowlist.size());
+
+  DataTable visdt;
+  visdt.AddDoubleColumn("mfc");
+  visdt.AddDoubleColumn("mtt");
+  visdt.AddIntegerColumn("jfreq");
+  visdt.AddIntegerColumn("numch");
+  visdt.AddFloatColumn("vre");
+  visdt.AddFloatColumn("vim");
+
+  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=0;
+  
+  char fname[1024];
+
+  TVector< uint_4 > chanum;
+  sprintf(fname, "%s/chanum.ppf",inoutpath.c_str());
+  {
+  PInPersist pic(fname);
+  pic >> chanum;
+  }
+
+  sa_size_t nrows = (imax-imin+1)/istep;\
+  sa_size_t kr=0;
+
+  for(int ifile=imin; ifile<=imax; ifile+=istep) {
+    sprintf(fname, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
+    cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
+    PInPersist pin(fname);
+    TMatrix< complex<r_4> > vismtx;
+    pin >> vismtx;
+
+    if (ifile==imin) {
+      sa_size_t ncols = vismtx.NCols();
+      cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles" 
+           << nrows << " NCols=NFreq=" << ncols << endl;
+      for(size_t j=0; j<rowlist.size(); j++)   vmtf[j].SetSize(nrows, ncols);
+    }
+    for(size_t j=0; j<rowlist.size(); j++) 
+      vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
+    kr++;
+     
+// Calcul moyenne dans des bandes en frequence 
+    FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
+  }
+
+  sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
+  cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;  
+  POutPersist po(fname);
+  char tagb[64];
+  for(size_t j=0; j<rowlist.size(); j++)  {
+    sprintf(tagb,"visft%d", chanum(rowlist[j]));
+    po << PPFNameTag(tagb) << vmtf[j];
+  }
+  cout << visdt;
+  sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
+  cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
+  POutPersist pod(fname);
+  pod << visdt;
+
+  return 0;
+}
+
+/* --Fonction-- */
+int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_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).real()/(r_4)(nmean_);
+        xnt_[5]=vismtx_(rv,jf).imag()/(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.;
+        for(sa_size_t jjf=jf; jjf<jf+djf_; jjf++) {
+          moyreal+=vismtx_(rv,jjf).real();
+          moyimag+=vismtx_(rv,jjf).imag();
+        }
+        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;
+}
+
+
+//-------------------------------------------------------------------- 
+//       Traitement fichiers produits par mcrd (V Jun09) 
+/* --Fonction-- */
+int DecodeProcVJun09(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;
+  if (narg>3) card=atoi(arg[3]); 
+  cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
+       << imin << "," << imax << "," << istep << " Card=" << card << endl;
+  cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
+  int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
+  return rc;
+}
+
+
+// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
+/* --Fonction-- */
 int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
 {
@@ -113 +309 @@
     else 
       sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
-    cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
+    cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
     PInPersist pin(fname);
     string tag1="specV1";
@@ -218 +414 @@
 
 
-
-
-/* --Fonction-- */
-int Usage(bool fgshort)
-{
-  cout << " --- svv2mtx.cc : Read PPF files produced by mcrd to make matrices BAORadio" << endl;
-  cout << " Usage:  mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1]" << endl; 
-  if (fgshort) {
-    cout << " mcrd -h for detailed instructions" << endl;
-    return 1;
-  }
-  cout << " InOutPath : Input/Output directory name " << endl;
-  cout << " Imin,Imax,IStep: Input PPF files sequence number \n"
-       << "    FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n" 
-       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
-       << " card=1 Ch12 , card=2 Ch34 " << endl;
-  return 1;
-}