Ajout fonctionalite remplissage cartes temps-frequence des Xcor, Reza 19/05/2015

visi2dtacx.cc

@@ -48,7 +48,8 @@ int Usage(bool fgea)
 {
   cout << " --- visi2dtacx.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
   if (fgea) cout << " visi2dtacx: Missing or wrong argument ! " << endl;
-  cout << " Usage: visi2dtacx InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
+  cout << " Usage: visi2dtacx [-tfm Nf] InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
+       << "  -tfm Nf : create time-frequncy maps of cross-correlation binning Nf frequencies \n" 
        << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
        << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
        << "  GainPPFFile: Input gains PPF file name \n"
@@ -72,26 +73,37 @@ int main(int narg, char* arg[])
 {
   if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
   if (narg<10) return Usage(true);
-  string path5 = arg[1];
-  string path6 = arg[2];
+  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<12) 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[3],"%d,%d",&Imin,&Imax);
+  sscanf(arg[offa+3],"%d,%d",&Imin,&Imax);
   Istep=1;
-  string gainfile=arg[4];
-  string outfile=arg[5];
+  string gainfile=arg[offa+4];
+  string outfile=arg[offa+5];
   // frequency range definition 
   double deltafreq=1.;
   double freq1=1375.;
   double freq2=1410.;
-  if ((narg>7)&&(strcmp(arg[6],"-")!=0))  freq1=atof(arg[6]);
-  if ((narg>8)&&(strcmp(arg[7],"-")!=0))  freq2=atof(arg[7]);
-  if ((narg>6)&&(strcmp(arg[8],"-")!=0))  deltafreq=atof(arg[8]);
+  int nargo=narg-offa;
+  if ((nargo>7)&&(strcmp(arg[offa+6],"-")!=0))  freq1=atof(arg[offa+6]);
+  if ((nargo>8)&&(strcmp(arg[offa+7],"-")!=0))  freq2=atof(arg[offa+7]);
+  if ((nargo>6)&&(strcmp(arg[offa+8],"-")!=0))  deltafreq=atof(arg[offa+8]);
   // time averaging interval definition, in term of visibility files
   int deltaIavg=10;
-  if ((narg>9)&&(strcmp(arg[9],"-")!=0))  deltaIavg=atoi(arg[9]);
+  if ((nargo>9)&&(strcmp(arg[offa+9],"-")!=0))  deltaIavg=atoi(arg[offa+9]);
   if (deltaIavg<1)  deltaIavg=1;
   int prtlev=1;
-  if (narg>10) prtlev=atoi(arg[10]);
+  if (nargo>10) prtlev=atoi(arg[offa+10]);
 
   cout << " visi2dtacx/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
        << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
@@ -99,6 +111,11 @@ int main(int narg, char* arg[])
     //       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
        <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
 
+  string TFMoutfile = "tfm_"+outfile;
+  if (FgTFM) {
+    cout << " visi2dtacx/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 
@@ -197,16 +214,28 @@ int main(int narg, char* arg[])
     double mttag;
     int cnt=0, cntnt=0, pcntnt=0;
     int I=0; 
-
+    //----- 6 TimeFrequency maps 
+    vector< TArray< complex<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;
+
     while (fgok) {
       fgok=vreader.ReadNext(vismtx, cfdate, mttag);
       if (fgok)  {
 	if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
 	  acsum.SetSize(8, vismtx.NCols());
 	  cxsum.SetSize(6, vismtx.NCols());
+	  if (FgTFM) {  //allocating time-frequency maps 
+	    tfmSX=(Imax-Imin)/Istep/deltaIavg;
+	    tfmSY=vismtx.NCols()/TFMfbin;
+	    cout << " visi2dtacx/Info: allocating Time-Freqncy maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
+	    for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) ); 
+	  }
 	}
 	if (I==0) {   // start filling a new DataTable row 
 	  dateobs=cfdate;
@@ -236,13 +265,23 @@ int main(int narg, char* arg[])
 	  }
 	  //---- On s'occupe des cross-correlations  1H-2H ... 3H-4H 
 	  TVector< complex<r_8> > vcx(acsum.NCols());
-	  for(int k=0; k<6; k++)  {
+	  for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
 	    for(sa_size_t c=0; c<vcx.Size(); c++) 
 	      vcx(c)=complex<r_8>((r_8)cxsum(k,c).real(), (r_8)cxsum(k,c).imag())*complex<r_8>(gaincx(k,c), 0.);
 	    cxdt[k] += vcx(freqs).Sum();    // integration en bande 1
 	    cxdt[k+6] += vcx(freqsB).Sum();    // integration en bande large B
 	    cxdt[k+12] += vcx(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
-	  }
+
+	    if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
+	      TArray< complex<r_4> > & tfmap = vtfm[k];
+	      for(sa_size_t jy=0; jy<tfmSY; jy++) {
+		complex<r_8> zz = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
+		tfmap(TFMtmidx, jy) = complex<r_4>(zz.real(), zz.imag());
+	      } 
+	    }  //----- end of time-frequency maps filling 
+	  }  //----- end of loop over the 6 Xcor 
+	  TFMtmidx++;
+
 	  // Moyennes dans la bande 1
 	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
 	  for(int k=0; k<8; k++)  acdt[k]*=fnorm; 
@@ -270,7 +309,8 @@ int main(int narg, char* arg[])
 	}
 	cnt++;
 	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
-	  cout << " visi2dtacx/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
+	  cout << " visi2dtacx/Info: fill-count="<<cntnt<<" fileCount="<<cnt<<" /Max="<<Imax
+	       <<" DateObs="<<dateobs<<endl;
 	  pcntnt=cntnt;
 	}
       } 
@@ -278,9 +318,20 @@ int main(int narg, char* arg[])
     cout << " visi2dtacx/Info: count="<<cnt<<" visimtx read "<<endl;
     dt.SetShowMinMaxFlag(true);
     cout << dt;
-    cout<<"  visi2dtacx/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
+    cout<<"  visi2dtacx/Info: Saving auto-corr=f(time) datatable to PPF file "<<outfile<<endl;
     POutPersist po(outfile);
     po<<dt; 
+
+    if (FgTFM) {   // --- renormalizing and saving time-frequency maps 
+      cout<<"  visi2dtacx/Info: Saving time-frequency maps to PPF file "<<TFMoutfile<<endl;
+      POutPersist potfm(TFMoutfile);
+      const char* tfm_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
+      for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
+	TArray< complex<r_4> > & tfmap = vtfm[k];
+	tfmap /= complex<r_4>((r_4)(1./(double)deltaIavg/(double)TFMfbin), 0.);
+	potfm << PPFNameTag(tfm_names[k]) << tfmap;
+      }
+    }
     //    resu.Update();
     cout << resu;   // Update est fait lors du print
   }