Changeset 3999 in Sophya for trunk/AddOn/TAcq/corrcrtadc.cc

Timestamp:

Jun 14, 2011, 10:21:16 PM (14 years ago)

Author:

ansari

Message:

amelioration prog calcul correlation ADC-CRT, Reza 14/06/2011

File:

• trunk/AddOn/TAcq/corrcrtadc.cc (modified) (11 diffs)

trunk/AddOn/TAcq/corrcrtadc.cc

@@ -4 +4 @@
 
 /* ---------------------------------------------------------- 
-   Projet BAORadio - (C) LAL/IRFU  2008-2010
+   Projet BAORadio - (C) LAL/IRFU  2008-2011
 
    Programme de lecture des fichiers DUMP ADC du CRT (Pittsburgh) 
    pour calcul de spectres / correlations 
    R. Ansari, C. Magneville   -  LAL/Irfu - Juin 2011
-   V : Mai 2009
    ---------------------------------------------------------- */
 
@@ -33 +32 @@
 #include "fftpserver.h"
 #include "utilarr.h"
+#include "histats.h"
 
 // include sophya mesure ressource CPU/memoire ...
@@ -43 +43 @@
 int Usage(bool fgshort=true);
 int DecodeProc(int narg, char* arg[]);
-int ProcessADCFiles(string& inoutpath, string& outname, int imin, int imax, int istep,  int jf1, int jf2, int nfreq);
+int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep,  int jf1, int jf2, int nfreq);
 int ComputeCorrel(TMatrix< complex<r_4> >& cfour,  TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel);
+int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep);
+
 //------------------------------------------------------------------------------------------------------------
 
@@ -50 +52 @@
 int Usage(bool fgshort)
 {
-  cout << " --- corrcrtadc.cc : Read CRT-ADC dump files to produce time-frequency\n"
-       << "     correlation matrices " << endl;
-  cout << " Usage: corrcrtadc InOutPath OutFile Imin,Imax,step [NumFreq1,NumFreq2,NBinFreq] \n" << endl;
+  cout << " --- corrcrtadc.cc : Read CRT-ADC dump files and process to produce correlation matrices \n"
+       << " Usage: corrcrtadc -corr/-hval InOutPath OutFile [AdcIds [Imin,Imax,step\n" << endl;
   if (fgshort) {
     cout << " corrcrtadc -h for detailed instructions" << endl;
     return 1;
   }
-  cout << " InOutPath : Input/Output directory name \n" 
+  cout << " -corr/-hval : Compute correlations or Fill histograms for ADC sample values \n" 
+       << " InOutPath : Input/Output directory name \n" 
+       << " ADCIds : ADC Id's to be processed (1,2,3,4)  \n" 
        << " OutFile : Output PPF file  name \n" 
        << " Imin,Imax,IStep: Input ADC dump files sequence number \n"
-       << "    FileNames=InOutPath/adcdumpNJFII.fits Imin<=II<=Imax II+=IStep , J=0,1,2 \n" 
-       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple" << endl;
+       << "    FileNames=InOutPath/adcdumpNJFII.fits Imin<=II<=Imax II+=IStep , J=ADCIds \n" << endl;
+  //       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple" << endl;
   return 1;
 }
@@ -101 +104 @@
 }
 
+
 //-------------------------------------------------------------------- 
 //       Traitement fichiers produits par vismfib  (V Nov09)
@@ -107 +111 @@
 {
   // Decodage des arguments et traitement 
-  string inoutpath = arg[0];  
-  string outname = arg[1];
+  string popt = arg[0];
+  bool fillhis=false;
+  if (popt=="-hval")  fillhis=true;
+  string inoutpath = arg[1];  
+  string outname = arg[2];
+  int aids[4]={0,1,-1,1};
+  if (narg>3) sscanf(arg[3],"%d,%d,%d,%d",aids,aids+1,aids+2,aids+3);
+  vector<int> adcids;
+  for(int ii=0;ii<4;ii++) 
+    if((aids[ii]>0)&&(aids[ii]<16))  adcids.push_back(aids[ii]);
   int imin=0;
   int imax=0;
   int istep=1;
-  sscanf(arg[2],"%d,%d,%d",&imin,&imax,&istep);
+  if (narg>4) sscanf(arg[4],"%d,%d,%d",&imin,&imax,&istep);
   int jf1=0;
   int jf2=0;
   int nfreq=0;
-  if (narg>3)
-    sscanf(arg[3],"%d,%d,%d",&jf1,&jf2,&nfreq);
-
-  cout << " ----- corrcrtadc/DecodeProc - Start - InOutPath= " << inoutpath << " OutFileName=" << outname << endl;
+  if (narg>5)
+    sscanf(arg[5],"%d,%d,%d",&jf1,&jf2,&nfreq);
+
+  cout << " ----- corrcrtadc/DecodeProc - Start " << ((fillhis)? " Fill ADC-val histogram" : " Compute Correlation") << endl;
+  cout <<  " InOutPath= " << inoutpath << " ADCIds: " ;
+  for(size_t ii=0; ii<adcids.size(); ii++) cout << adcids[ii] << " , ";
+  cout << " (NbADC=" << adcids.size() << ") OutFileName=" << outname << endl;
   cout << " IMin,Max,Step="    << imin << "," << imax << "," << istep 
        << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
-  int rc=ProcessADCFiles(inoutpath, outname, imin, imax, istep, jf1, jf2, nfreq);
+  int rc=0;
+  if (fillhis) rc = ADCFiles2Histo(inoutpath, adcids, outname, imin, imax, istep);
+  else rc=ProcessADCFiles(inoutpath, adcids, outname, imin, imax, istep, jf1, jf2, nfreq);
+  
   return rc;
 }
 
-// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
-/* --Fonction-- */
-int ProcessADCFiles(string& inoutpath, string& outname, int imin, int imax, int istep, int jf1, int jf2, int nfreq) 
-{
-  Timer tm("ProcessADCFiles");  
-
-  // Taille des fichiers 4 voies = 1024*1024*128 ===> 32 x 1024 x 1024 samples / voie
-  // On lit par paquet de 4096 = 4 x 1024 / voie --> ADCRDBLKSZ=16x1024=16384
+// Taille des fichiers 4 voies = 1024*1024*128 ===> 32 x 1024 x 1024 samples / voie
+// On lit par paquet de 4096 = 4 x 1024 / voie --> ADCRDBLKSZ=16x1024=16384
 #define ADCFLEN 4096 
 #define ADCRDBLKSZ 16384
 #define ADCNREAD 8192 
-#define NBADC 2
-#define NFREQUSED 2000
+
+static int PRTMODULO = 250;
+
+// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
+/* --Fonction-- */
+int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep, int jf1, int jf2, int nfreq) 
+{
+  Timer tm("ProcessADCFiles");  
+
+#define NFREQUSED 2040
+#define NFREQUSESTART 5
+#define NFREBIN 4 
+
+  sa_size_t NFREQREBIN = NFREQUSED/NFREBIN;
+  sa_size_t NBADC = adcids.size();
 
   sa_size_t  NCHAN = 4*NBADC;
@@ -145 +170 @@
   TMatrix< complex<r_8> > correl(NCORREL,NFREQUSED);
   TMatrix< complex<r_8> > autocorrel(NCHAN,NFREQUSED);
-
-  FILE* fip[NBADC]={NULL,NULL};
+  double nsumcor=0.;
+
+  vector<FILE*> fip; 
+  vector<int> adcchans;
+  for(size_t ka=0; ka<NBADC; ka++)  { 
+    fip.push_back(NULL);
+    adcchans.push_back(adcids[ka]*4+ka);
+  }
+
   FFTPackServer ffts(false);
   
@@ -155 +187 @@
 
   signed char rdbuff[ADCRDBLKSZ]; 
-  int prtmodulo=50;
   for(int ifile=imin; ifile<=imax; ifile+=istep) {
     for(int ka=0; ka<NBADC; ka++) {
-      sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),ka,ifile);
+      sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
       fip[ka]=fopen(fname,"rb");
-      if (fip[ka]==NULL) cout << "ProcessADCFiles/ERROR opening file " << fname << endl;
-      else  cout << "ProcessADCFiles/Info - file " << fname << " opened ..." << endl;
+      if (fip[ka]==NULL) 
+        cout << "ProcessADCFiles/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
+      else  
+        cout << "ProcessADCFiles/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
     }
 
@@ -171 +204 @@
           for(sa_size_t ls=0; ls<ADCFLEN; ls++) sigadc(ls)=(float)rdbuff[ls*4+jac];
           ffts.FFTForward(sigadc, foursig);
-          //      cfour.Row(irc)=foursig(Range(400,1679)).Transpose();  irc++;
-          cfour.Row(irc)=foursig(Range(10,2009)).Transpose();  irc++;
-
+          cfour.Row(irc)=foursig(Range(NFREQUSESTART,NFREQUSESTART+NFREQUSED-1)).Transpose();  irc++;
         }
       }
       ComputeCorrel(cfour, correl, autocorrel);
-      if (m%prtmodulo==0)  cout << " ProcessADCFiles/Info Done read+FFT+correl m=" << m << " / Max=" << ADCNREAD << endl;
+      nsumcor++;
+      if (m%PRTMODULO==0)  cout << " ProcessADCFiles/Info Done read+FFT+correl m=" << m << " / Max=" << ADCNREAD << endl;
     }
     for(int ka=0; ka<NBADC; ka++)  fclose(fip[ka]);
   }
 
+
+  correl /= nsumcor;
+  autocorrel /= nsumcor;
+
+  TMatrix< complex<r_8> > correbin(NCORREL,NFREQREBIN);
+  sa_size_t icc=0;
+  for(sa_size_t i=0; i<correbin.NCols(); i++) {
+    for(sa_size_t j=0; j<NFREBIN; j++) { 
+      correbin.Column(i) += correl.Column(icc); 
+      icc++;
+    }
+  }
 
   sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
@@ -188 +232 @@
   po << PPFNameTag("correlmtx") << correl; 
   po << PPFNameTag("acorrmtx") << autocorrel; 
+  po << PPFNameTag("rebincormtx") << correbin; 
 
   cout << "ProcessADCFiles:  finished FFT + correlation computing " << endl;
+
+  cout << " --------------- Channel/Correlation number association --------------- " << endl;
+  cout << " CorreNumber- axb  (MtxRow) :  (Ca, Cb)  -- ADCChan(c,d) " << endl; 
+  cout << " ----------------------------------------------------------------------- " << endl;
+  sa_size_t jcr=0;
+  for(sa_size_t j1=0; j1<cfour.NRows(); j1++) {
+    for(sa_size_t j2=j1; j2<cfour.NRows(); j2++) {
+      cout <<  jcr << " : (" << j1 << "," << j2 << ") -> ADCChans (" << adcchans[j1] << "," << adcchans[j2] << ")" << endl;
+      jcr++;
+    }
+  }
+  cout << " ----------------------------------------------------------------------- " << endl;
+
   return 0;
 }
@@ -208 +266 @@
 }
 
+/* --Fonction-- */
+int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep)
+{
+  Timer tm("ADCFiles2Histo");  
+
+  sa_size_t NBADC = adcids.size();
+
+  sa_size_t  NCHAN = 4*NBADC;
+  vector<Histo *> vhist;
+  for(int ka=0; ka<NCHAN; ka++) 
+    vhist.push_back(new Histo(-128.5,128.5,257) );
+
+  vector<FILE*> fip; 
+  for(size_t ka=0; ka<NBADC; ka++)  fip.push_back(NULL);
+    
+  char fname[512];
+
+  signed char rdbuff[ADCRDBLKSZ]; 
+  for(int ifile=imin; ifile<=imax; ifile+=istep) {
+    for(int ka=0; ka<NBADC; ka++) {
+      sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
+      fip[ka]=fopen(fname,"rb");
+      if (fip[ka]==NULL) 
+        cout << "ADCFiles2Histo/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
+      else  
+        cout << "ADCFiles2Histo/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
+    }
+
+    for(int m=0; m<ADCNREAD; m++) {
+      sa_size_t irc=0;
+      for(int ka=0; ka<NBADC; ka++) {
+        fread(rdbuff,1,ADCRDBLKSZ,fip[ka]);
+        for(int jac=0; jac<4; jac++) {
+          for(sa_size_t ls=0; ls<ADCFLEN; ls++) vhist[irc]->Add((r_8)rdbuff[ls*4+jac]);
+          irc++;
+        }
+      }
+      if (m%PRTMODULO==0)  cout << " ADCFiles2Histo/Info Done FillHist m=" << m << " / Max=" << ADCNREAD << endl;
+    }
+    for(int ka=0; ka<NBADC; ka++)  fclose(fip[ka]);
+  }
+
+  sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
+  cout << "ADCFiles2Histo: Opening file " << fname << " for writing ADC value histograms " << endl;  
+  POutPersist po(fname);
+  char nametag[64];
+  for(size_t ka=0; ka<NCHAN; ka++) {
+    sprintf(nametag,"adcvalC%d",ka);
+    po << PPFNameTag("nametag") << *(vhist[ka]); 
+  }
+  for(size_t ka=0; ka<NCHAN; ka++) delete vhist[ka];
+
+  cout << "ADCFiles2Histo:  finished ADC value histogram fill " << endl;
+  return 0; 
+}