Decodage -reorderfreq ds P4AnaParams , Reza 28/2/2017

README

 -------------------------------------------------------------------
   PAON-4 visibility analysis programs 
-       R. Ansari, J.E. Campagne, C. Magneville   -  Feb 2015 
+       R. Ansari, J.E. Campagne, C. Magneville   -  Feb 2015
+                                                    Updated Feb 2017 
  -------------------------------------------------------------------
 
 
@@ -17,16 +18,25 @@ List of files:
     Updated April 2015: computes also a frequency dependent gain gain(nu) using a gain 
     function developed by Christophe (cmv) for cluster data analysis 
     Updated July 2016: use -m option to read the mean FFT coef matrix
-- visi2ntac.cc : a simple main program to read visibilitiy files and produce a NTuple 
-  with total power as a function of time for the 8 PAON4 auto-correlations 
-    Updated April 2015: reads the gain g(nu) file produced by rdvisip4.cc , corrects the 
-    auto-correlation signals before integrating in the specified bands 
+
+- p4autils.h p4autils.cc : utility classes and functions, in particular argument decoder class
+    for different programs 
+    
 - visi2dtacx.cc : a main program to read visibilitiy files and produce a datatable 
   with auto-correlation and 6 cross correlations (1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H) 
   as a function of time for three frequency bins - Can also make a time-frequency map of 
   the 6 cross-correlations. The visibilities are corrected for the electronic gain, using the 
   gain curve g(nu) computed by rdvisip4.cc 
+
 - visi2tmfreq.cc : a main program to read visibilitiy files and produce time-frequency maps 
-  for the 8 auto-correlations, corrected for the gain, using the gain curve g(nu) computed by rdvisip4.cc 
+  for the 8 auto-correlations, corrected for the gain, using the gain curve g(nu) computed by rdvisip4.cc
+  Updated Feb 2017 : computes and fills an NTuple with mean/sigma for the 8 auto-correlations signals
+  
 - msvis2dt.cc (sep 2015) : a program to read meanspecII.ppf (firmware RAW, monitoring output) 
 and vismtxII.ppf files ( 4 feeds) to create datatable and time-frequency maps  
+
+- visi2ntac.cc :  DEPRECATED   - replaced by visi2tmfreq
+  a simple main program to read visibilitiy files and produce a NTuple 
+  with total power as a function of time for the 8 PAON4 auto-correlations 
+    Updated April 2015: reads the gain g(nu) file produced by rdvisip4.cc , corrects the 
+    auto-correlation signals before integrating in the specified bands 

makefile

@@ -11,7 +11,7 @@ all : rdvisip4  visi2ntac  visi2dtacx  visi2tmfreq p4conv2fits msvis2dt
 clean :
 	rm -f $(EXE)/rdvisip4  $(EXE)/visi2ntac  $(EXE)/visi2dtacx  $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt
 	rm -f $(OBJ)/rdvisip4.o  $(OBJ)/visi2ntac.o  $(OBJ)/visi2dtacx.o  $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o
-	rm -f $(OBJ)/visip4reader.o
+	rm -f $(OBJ)/p4autils.o  $(OBJ)/visip4reader.o
 
 ######
 ## programme de remplissage de DataTable avec cross-correlation fct du temps
@@ -40,10 +40,10 @@ $(OBJ)/visi2ntac.o : visi2ntac.cc visip4reader.h
 visi2tmfreq : $(EXE)/visi2tmfreq
 	echo '---visi2tmfreq made'
 
-$(EXE)/visi2tmfreq : $(OBJ)/visi2tmfreq.o $(OBJ)/visip4reader.o 
-	$(CXXLINK) -o $(EXE)/visi2tmfreq $(OBJ)/visi2tmfreq.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST) 
+$(EXE)/visi2tmfreq : $(OBJ)/visi2tmfreq.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o 
+	$(CXXLINK) -o $(EXE)/visi2tmfreq $(OBJ)/visi2tmfreq.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST) 
 
-$(OBJ)/visi2tmfreq.o : visi2tmfreq.cc visip4reader.h
+$(OBJ)/visi2tmfreq.o : visi2tmfreq.cc visip4reader.h p4autils.h
 	$(CXXCOMPILE)  -o $(OBJ)/visi2tmfreq.o visi2tmfreq.cc
 
 ######
@@ -80,7 +80,10 @@ $(EXE)/msvis2dt : $(OBJ)/msvis2dt.o $(OBJ)/msvis2dt.o
 $(OBJ)/msvis2dt.o : msvis2dt.cc 
 	$(CXXCOMPILE)  -o $(OBJ)/msvis2dt.o msvis2dt.cc
 
+######
+$(OBJ)/p4autils.o : p4autils.cc p4autils.h
+	$(CXXCOMPILE)  -o $(OBJ)/p4autils.o p4autils.cc
+
 ######
 $(OBJ)/visip4reader.o : visip4reader.cc visip4reader.h
 	$(CXXCOMPILE)  -o $(OBJ)/visip4reader.o visip4reader.cc
-

p4autils.cc

@@ -58,14 +58,17 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
       if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
       gain_gnu_file_=arg[2];   arg+=2; narg-=2; 
     }
-    else if (fbo=="-o") {
-      if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
-      outfile_=arg[2];   arg+=2; narg-=2; 
-    }
     else if (fbo=="-nmean") {
       if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
       Nmean_=atoi(arg[2]);   arg+=2; narg-=2; 
     }
+    else if (fbo=="-reorderfreq") {
+      fgreorderfreq_=true;   arg++; narg--; 
+    }
+    else if (fbo=="-o") {
+      if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
+      outfile_=arg[2];   arg+=2; narg-=2; 
+    }
     else if (fbo=="-tfm") {
       if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
       int ja,jb;
@@ -87,17 +90,18 @@ int P4AnaParams::UsageOptions()
 {
   cout<< " P4AnaParams::UsageOptions(): -inrange -in5 -in6 -gnu -o -nmean -tfm -freq -prt \n"
       << "    [-inrange Imin,Imax[,Istep] ]  [-in5 InputPath_BAO5]   [-in6 InputPath_BAO6]  \n"
-      << "    [-o outfilename] [-gnu gain_filename] [-nmean N]  [-prt level] \n"
-      << "    [-tfm timebin,freqbin]  [-freq f0,df] \n"
+      << "    [-o outfilename] [-gnu gain_filename] [-nmean N] [-reorderfreq] \n"
+      << "     [-tfm timebin,freqbin]  [-freq f0,df] [-prt level] \n"
       << " -inrange Imin,Imax[,Istep]  define range of visibilite matrices to be processed (Imin<=I<=Imax with step Istep) \n"
       << " -in5/6 InputPath_BAO5/6 : input file path for BAO5 and BAO6 \n"
       << " -o outfilename : specify output file name \n"
       << " -gnu gain_filename: gain g(nu) file name \n"
       << " -nmean N: compute average (and sigma...) using N visibility matrices \n"
-      << " -prt level: define print level  \n"
+      << " -reorderfreq: activate frequency reordering in VisiP4Reader, suitable for FFT firmware \n"
       << " -tfm timebin,freqbin : time-frequency maps, averaging timebin=nb of visi matrices, freqbin frequency components \n"
       << " -freq f0,df : definition of a frequency band centered on f0, with total width df (in MHz) \n"
       << "     several frequency bands can be specified by mutiple -freq options \n"
+      << " -prt level: define print level  \n"
       << endl;
   return 0;
 }

rdvisip4.cc

@@ -89,7 +89,7 @@ int main(int narg, char* arg[])
 
   string path5 = arg[1+off];
   string path6 = arg[2+off];
-  int Imin,Imax,Istep; 
+  int Imin=0,Imax=10,Istep=1; 
   sscanf(arg[3+off],"%d,%d,%d",&Imin,&Imax,&Istep);
   string outfile=arg[4+off];
   string gainfile=arg[5+off];

visi2tmfreq.cc

@@ -8,7 +8,7 @@
    visi2tmfreq.cc: programme de lecture des fichiers matrices de 
    visibilites de PAON4 produits mfacq, creation d'une tableau 
    temps-frequence pour chaque auto-correlation 
-   R. Ansari, T. Etourneau   -  LAL
+   R. Ansari  -  LAL   (stage de T. Etourneau)
    V1 : Juin 2015 
    ---------------------------------------------------------- */
 
@@ -38,77 +38,53 @@
 #include "timing.h"
 
 // include lecteur de fichiers visibilites 
+#include "p4autils.h"
 #include "visip4reader.h"
 
  
-//--------------------------- Fonctions de ce fichier   ------------------- 
-int Usage(bool fgshort=true);
-// int DecodeArgs(int narg, char* arg[]);
+//--------------------------- Fonctions / classe de ce fichier   ------------------- 
+NTuple CreateAC_MeanSigNTuple(P4AnaParams const& p);
 
 /* --Fonction-- */
-int Usage(bool fgea)
+int Usage()
 {
-  cout << " --- visi2tmfreq.cc : Read PPF files produced by mfacq to create time-frequency maps\n" << endl;
-  if (fgea) cout << " visi2tmfreq: Missing or wrong argument ! " << endl;
-  cout << " Usage: visi2tmfreq InPathBAO5 InPathBAO6 Imin,Imax  GainPPFFile OutPPFFile [DeltaIAvg=60] [NBinFreq=4] [PrintLev=0] \n"
-       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
-       << "  Imin,Imax : read files for Imin<=iii<=Imax  \n"
-       << "  GainPPFFile: Input gains PPF file name \n"
-       << "  OutPPFFile: Output PPF file name which will contain the 8 auto-correlation time-frequency arrays \n"
-       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=60 \n"
-       << "  NBinFreq: The binning in frequency in the computed time-frequency maps (def=4 -> 250 kHz) \n"
-       << endl;
-  /*
-  if (fgshort) {
-    cout << " rdvisip4 -h for detailed instructions" << endl;
-    return 1;
-  }
-  */
+  cout << " --- visi2tmfreq.cc : Read PPF files produced by mfacq to create time-frequency maps and autocorre mean/sigma NTuple \n" << endl;
+  cout << " Usage: visi2tmfreq [-arguments] [OutPPFFile] \n"<<endl;
+  P4AnaParams::UsageOptions();
+  cout<< endl;
   return 1;
 }
 
 //----------------------------------------------------
 //----------------------------------------------------
-int main(int narg, char* arg[])
+int main(int narg, const char* arg[])
 {
-  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
-  if (narg<6) return Usage(true);
-  string path5 = arg[1];
-  string path6 = arg[2];
-  int Imin,Imax,Istep=1; 
-  sscanf(arg[3],"%d,%d",&Imin,&Imax);
-  string gainfile=arg[4];
-  string outfile=arg[5];
-  Istep=1;
-
-  // taille de fenetre de moyenne en temps
-  int deltaIavg=60;
-  if ((narg>6)&&(strcmp(arg[6],"-")!=0))  deltaIavg=atoi(arg[6]);
-  // taille de fenetre de moyenne en frequence (en bin de frequence, cad 61 kHz) 
-  sa_size_t TFMfbin=4;
-  if ((narg>7)&&(strcmp(arg[7],"-")!=0))  TFMfbin=atoi(arg[7]);
-  int prtlev=0;
-  if (narg>8) prtlev=atoi(arg[8]);
-
-  cout << " visi2tmfreq/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
-       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<'\n'
-       << " GainFile: "<<gainfile<<" OutFile: "<<outfile
-       << " DeltaIAvg="<<deltaIavg<<" NBinFreq="<<TFMfbin<<" PrtLev="<<prtlev<<endl;
+  if ((narg<8)||((narg>1)&&(strcmp(arg[1],"-h")==0)) )  return Usage();
   HiStatsInitiator _inia;
   //   TArrayInitiator  _inia;
 
   int rc = 0;
   try {
+    // Decoding parameters 
+    P4AnaParams params;
+    params.DecodeArgs(narg, arg);
+    string outfile=params.outfile_;
+    if (params.lastargs_.size()>0)  outfile=params.lastargs_[0];
+    
+    cout << " visi2tmfreq/Info: \n"<<params<<endl;
+    int prtlev=params.prtlev_;
+    
     ResourceUsage resu;
     //  Numero des lignes des auto-correlations dans la matrice des visibilites 
     sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
  
-    cout << " visi2tmfreq/Info: reading input gain file"<<gainfile<<" ..."; 
+    cout << " visi2tmfreq/Info: reading input gain file"<<params.gain_gnu_file_<<" ...";
+    /* ---- NOTE/Remarque: mettre la lecture des gains dans une fonction / classe --- */
     Matrix gains;
     Vector gn;
     double GAINTHR = 0.005;  // limite inferieur du gain avant faire 1/G
     {
-      PInPersist pin(gainfile);
+      PInPersist pin(params.gain_gnu_file_);
       pin>>PPFNameTag("gains")>>gains;
       pin>>PPFNameTag("gn")>>gn;
       for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
@@ -122,74 +98,143 @@ int main(int narg, char* arg[])
 
     // Creation de l'instance de la classe lecteur des fichiers visimtxXXX.ppf 
     vector<string> paths; 
-    paths.push_back(path5); 
-    paths.push_back(path6);
-    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
-    vreader.setPrintLevel(prtlev);
+    paths.push_back(params.inpath5_); 
+    paths.push_back(params.inpath6_);
+    VisiP4Reader vreader(paths, params.Imin_, params.Imax_, params.Istep_, params.fgreorderfreq_);
+    vreader.setPrintLevel(params.prtlev_);
 
     bool fgok=true;
     TMatrix< complex<r_4> > vismtx;
     TMatrix< complex<r_4> > acsum;
+    TMatrix< r_4 > gcac;   // gain corrected Auto-corelation signals 
 
     double acdt[32];   // les 4*8=32 valeurs d'autocorrelation pour remplissage dans la table 
     complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 
 
-    TimeStamp dateobs, cfdate;
-    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
+    TimeStamp dateobs, cfdate, firstdate;
     double mttag;
     int cnt=0, pcnttm=0;
-    int I=0; 
+    // Index pour le remplissage de la carte temps-frequence 
+    int ITFM=0;
+    int TFM_DeltaI=params.TFMtimebin_;
     //----- 8 TimeFrequency maps of auto-correlation 
     vector< TArray< r_4 > > vtfm;
-    
+
+
     //---- for the time-freqency map filling    
     sa_size_t TFMtmidx=0;
     sa_size_t tfmSX, tfmSY;
 
+    /*---- NOTE/REMARQUE (reza+cmv): Regrouper le creation/calcul/remplissage du NTuple dans une classe 
+               ----  creation/calcul/remplissage NTuple - Partie I                ---- */
+    // Index pour le remplissage du NTuple 
+    int JNT=0, JJNT=0;
+    int JNT_DeltaJ=params.Nmean_;
+    NTuple nt=CreateAC_MeanSigNTuple(params);
+    double xnt[100];
+    double sump[50];  //  will contain sum of the power for each auto-correlation and each frequency band
+    double sumsqp[50];  //  will contain sum of the square of power  (p^2) for each auto-correlation and each frequency band 
+    int nbands=params.fbands_.size();
+    if (nbands>5) nbands=5;
+    bool fgnt=false;
+    if (nt.NRows()>0)  fgnt=true;  // il faut remplir le NTuple.
+    /*   ---- FIN  creation/calcul/remplissage NTuple - Partie I     ---- */
+
+    
     while (fgok) {  // boucle de lecture des fichiers d'entree 
       fgok=vreader.ReadNext(vismtx, cfdate, mttag);
       if (!fgok) continue; 
-	if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
-	  acsum.SetSize(8, vismtx.NCols());
-	  // Allocating time-frequency maps 
-	  tfmSX=vismtx.NCols()/TFMfbin;   // frequence selon l'axe X 
-	  tfmSY=(Imax-Imin)/Istep/deltaIavg;  // temps selon l'axe Y 
-	    cout << " visi2tmfreq/Info: allocating Time-Freqncy maps : Freq->NX= "<<tfmSX<<" x Time->NY="<<tfmSY<<endl;
-	    for(int k=0; k<8; k++) vtfm.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
-	  }
-	if (I==0) {   // start filling a new DataTable row 
-	  dateobs=cfdate;
-	  if (prtlev>0) 
-	    cout << "visi2tmfreq/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
-	  acsum = complex<r_4>(0.,0.);
-	}
-	//   sum (integration) along the time axis 
-	for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
+      if (cnt==0)  {
+	firstdate=cfdate;
+	//resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
+	acsum.SetSize(8, vismtx.NCols());
+	gcac.SetSize(8, vismtx.NCols());
+	// Allocating time-frequency maps 
+	tfmSX=vismtx.NCols()/params.TFMfreqbin_;   // frequence selon l'axe X 
+	tfmSY=(params.Imax_-params.Imin_)/params.Istep_/TFM_DeltaI;  // temps selon l'axe Y 
+	cout << " visi2tmfreq/Info: allocating Time-Freqncy maps : Freq->NX= "<<tfmSX<<" x Time->NY="<<tfmSY<<endl;
+	for(int k=0; k<8; k++) vtfm.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
+      }
+      if (ITFM==0) {   // start filling a ne time bin in the time-frequency map 
+	dateobs=cfdate;
+	if (prtlev>0) 
+	  cout << "visi2tmfreq/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
+	acsum = complex<r_4>(0.,0.);
+      }
+      //----- extracting the autocorrelation information 
+      for(int k=0; k<8; k++)  {   // boucle sur les auto-correlations
+	acsum.Row(k) += vismtx.Row(KVAC[k]);     //  sum (integration) along the time axis  
+	/// gcac.Row(k) = ComplexMathArray< r_4>::real(vismtx.Row(KVAC[k]), true);
+	for(sa_size_t c=0; c<gcac.NCols(); c++)  gcac(k,c)=vismtx(KVAC[k], c).real()*gains(k,c);
+      }
 
-      //   sum (integration) along the time axis 
-      for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
-      I++;    // incrementing DeltaTime counter 
-      if (I==deltaIavg) {
+      //------------- calcul/rempilssage de de la matrice temps-frequence 
+      ITFM++;    // incrementing DeltaTime counter 
+      if (ITFM==TFM_DeltaI) {
 	Vector vac(acsum.NCols());
 	// facteur de normalisation pour que notre carte temps-frequence corresponde a la moyenne et pas la somme 
-	double facnorm=1./(double)(TFMfbin*deltaIavg);  
+	double facnorm=1./(double)(params.TFMfreqbin_*TFM_DeltaI);  
 	for(int k=0; k<8; k++)  {  // loop over 8 auto-corr
 	  // correct for gain(nu)
 	  for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
 	  //filling time-frequency maps  
 	  TArray< r_4 > & tfmap = vtfm[k];
 	  for(sa_size_t ix=0; ix<tfmSX; ix++) 
-	    tfmap(ix, TFMtmidx) = (vac( Range(ix*TFMfbin, (ix+1)*TFMfbin-1) ).Sum())*facnorm;
+	    tfmap(ix, TFMtmidx) = (vac( Range(ix*params.TFMfreqbin_, (ix+1)*params.TFMfreqbin_-1) ).Sum())*facnorm;
 	}  //---- end of loop over 8 auto-corr
 	TFMtmidx++;
-	I=0; 
+	ITFM=0; 
       }
       cnt++;
       if ((TFMtmidx>0)&&(TFMtmidx%10==0)&&(TFMtmidx>pcnttm)) {
-	cout << " visi2tmfreq/Info: TimeIndex="<<TFMtmidx<<" fileCount="<<cnt<<" /Max="<<Imax
+	cout << " visi2tmfreq/Info: TimeIndex="<<TFMtmidx<<" fileCount="<<cnt<<" /Max="<<params.Imax_
 	     <<" DateObs="<<dateobs<<endl;
 	pcnttm=TFMtmidx;
       }
+      if (! fgnt)  continue;  // pas de NTuple a remplir
+    /*---- NOTE/REMARQUE (reza+cmv): Regrouper le creation/calcul/remplissage du NTuple dans une classe 
+               ----  creation/calcul/remplissage NTuple - Partie II                ---- */
+      // pour remplissage NTuple puissance moyenne et sigma 
+      if (JNT==0) {
+	dateobs=cfdate;
+	for(int j=0; j<nbands; j++) {   // boucle sur les bandes en frequence 
+	  for(int k=0; k<8; k++) {    // boucle sur les 8 voies d'auto-correlation
+	    sump[j*8+k]=sumsqp[j*8+k]=0.;
+	  }
+	}
+      }
+      for(int j=0; j<nbands; j++) {   // boucle sur les bandes en frequence
+	sa_size_t JFmin=params.fbands_[j].jfmin_;
+	sa_size_t JFmax=params.fbands_[j].jfmax_;
+	Range freqs(JFmin, JFmax);
+	double fnorm=1./(double)(JFmax-JFmin+1);
+	for(int k=0; k<8; k++) {    // boucle sur les 8 voies d'auto-correlation
+	  // puissance moyenne dans la bande
+	  TVector< r_4 > acors=gcac.Row(k);
+	  double signal_in_band=(acors.SubVector(freqs).Sum())*fnorm;
+	  sump[j*8+k]+=signal_in_band;
+	  sumsqp[j*8+k]+= (signal_in_band*signal_in_band);
+	}
+      }
+      JNT++;
+      if (JNT==JNT_DeltaJ)  {  // remplissage du NTuple 
+	double duration=TimeStamp::TimeDifferenceSeconds(cfdate,dateobs);
+	xnt[0]=JJNT;   xnt[1]=dateobs.DaysPart()-firstdate.DaysPart();
+	xnt[2]=dateobs.SecondsPart()+0.5*duration;   xnt[3]=TimeStamp::TimeDifferenceSeconds(dateobs,firstdate)+0.5*duration;
+	double fnorm=(1./(double)JNT_DeltaJ);
+	int kxnt=4;
+	for(int j=0; j<nbands; j++) {   // boucle sur les bandes en frequence
+	  for(int k=0; k<8; k++) {    // boucle sur les 8 voies d'auto-correlation
+	    sump[j*8+k] *= fnorm;
+	    xnt[kxnt]=sump[j*8+k];  kxnt++;
+	    sumsqp[j*8+k] *= fnorm;   sumsqp[j*8+k] -= (sump[j*8+k]*sump[j*8+k]);
+	    xnt[kxnt]=((sumsqp[j*8+k]>=0.)?sqrt(sumsqp[j*8+k]):-1.);  kxnt++;
+	  }
+	}
+	nt.Fill(xnt);
+	JNT=0;
+      }
+    /*   ---- FIN  creation/calcul/remplissage NTuple - Partie II     ---- */
 
     }    //--- fin de boucle de lecture des fichiers d'entree 
 
@@ -201,6 +246,11 @@ int main(int narg, char* arg[])
       TArray< r_4 > & tfmap = vtfm[k];
       potfm << PPFNameTag(tfm_names[k]) << tfmap;
     }
+    if (fgnt) {
+      cout<<" visi2tmfreq/Info: saving AutoCorrelation mean/sigma NTuple as NTAC to PPF file "<<outfile<<endl;
+      potfm << PPFNameTag("NTAC") << nt;
+      cout<<nt; 
+    }
     //    resu.Update();
     cout << resu;   // Update est fait lors du print
   }
@@ -224,3 +274,43 @@ int main(int narg, char* arg[])
 
 }
 
+
+/*---- NOTE/REMARQUE (reza+cmv): Regrouper le creation/calcul/remplissage du NTuple dans une classe 
+    ----  creation/calcul/remplissage NTuple - Partie III                ---- */
+
+/*-- Fonction --*/
+NTuple CreateAC_MeanSigNTuple(P4AnaParams const& p)
+{
+  if (!p.fgdtable_ || p.fbands_.size()<1)  {
+    cout<<" CreateAC_MeanSigNTuple() no frequency band specified - NTuple not created..."<<endl;
+    NTuple ntv;
+    return ntv;
+  }
+  int nbands=p.fbands_.size();
+  if (p.fbands_.size()>5) {
+    cout<<" CreateAC_MeanSigNTuple()  maximum 5 frequency bands can be filled in the NTuple..."<<endl;
+    nbands=5;
+  }
+  int cnt=0;
+  int L=16;  // longueur des noms
+  char buff[1600];
+  char* ntnames[100];
+  ntnames[cnt]=buff+cnt*L;  strcpy(ntnames[cnt], "k");  cnt++;
+  ntnames[cnt]=buff+cnt*L;  strcpy(ntnames[cnt], "day");  cnt++;
+  ntnames[cnt]=buff+cnt*L;  strcpy(ntnames[cnt], "hour");  cnt++;
+  ntnames[cnt]=buff+cnt*L;  strcpy(ntnames[cnt], "timeoff");  cnt++;
+  DVList info;
+  char vname[32];  
+  for(int j=0; j<nbands; j++) {   // boucle sur les bandes en frequence
+    sprintf(vname,"Freq%d",j+1);   info[vname]=p.fbands_[j].f0_;
+    sprintf(vname,"DeltaFreq%d",j+1);   info[vname]=p.fbands_[j].df_; 
+    for(int k=0; k<8; k++) {    // boucle sur les 8 voies d'auto-correlation
+      ntnames[cnt]=buff+cnt*L;  sprintf(ntnames[cnt],"p%df%d",k,j+1);   cnt++;
+      ntnames[cnt]=buff+cnt*L;  sprintf(ntnames[cnt],"sigp%df%d",k,j+1);   cnt++;
+    }
+  }
+  NTuple ntt(cnt, ntnames, 384, false);    // float numbers have enough precision for us
+  ntt.Info() = info;
+  return ntt;
+}
+/*   ---- FIN  creation/calcul/remplissage NTuple - Partie III     ---- */