Changeset 3943 in Sophya for trunk/AddOn/TAcq/brproc.cc

Timestamp:

Feb 1, 2011, 9:30:20 AM (15 years ago)

Author:

campagne

Message:

release.txt

File:

• trunk/AddOn/TAcq/brproc.cc (modified) (16 diffs)

trunk/AddOn/TAcq/brproc.cc

@@ -43 +43 @@
     fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan), 
     nbpaq4mean_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()), 
-    nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()) 
+    nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres())
 {
   setNameId("meanSpecCalc",1);
-  BRPaquet paq(memgr_.PaqSize());
+
+  uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
+  const char* venvp = NULL; 
+  venvp=getenv("BRANA_NBYTECUT");
+  if (venvp!=NULL){
+    nb_octets_entrop = atoi(venvp);
+    cout << "BRMeanSpecCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
+  }
+
+  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
+
   if (fgsinglechannel_) {
     mspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); 
@@ -70 +80 @@
   }
 
+  //
+  
+
   // Definition des tailles de fenetres de spectres, etc ...
   SetSpectraWindowSize();
-  SetMaxNbSepcWinFiles();
+  SetMaxNbSpecWinFiles();
   nbtot_specwin_=0;
   SetVarianceLimits();
+  SetNumberOfBands();
 
   ofsdtp_=NULL;
   dtp_=NULL;
-  xnt_=new double[nchan*2];  // CHECK : ATTENTION la taille depend de nombre de colonne du NTuple !
+  
+  //  cout << "(JEC) creation tuple de " << nchan*(1+numberOfBands_) << " doubles " << endl;
+  //  xnt_=new double[nchan*(1+numberOfBands_)];  // CHECK : ATTENTION la taille depend de nombre de colonne du NTuple !
+  xnt_=NULL;
 }
 
@@ -84 +101 @@
 BRMeanSpecCalculator::~BRMeanSpecCalculator()
 {
+  cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
   uint_8 npqm=0;
   for(size_t i=0; i<nbpaq4mean_.size(); i++)  npqm+=nbpaq4mean_[i];
   if (npqm>nmean_*nbpaq4mean_.size()/10)  SaveMeanSpectra();
-  cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
   for(size_t i=0; i<nbadpaq_.size(); i++) {
     cout << " Channel " << i << "  NBadPaq=" << nbadpaq_[i] << " / TotNbPaq=" << totnbpaq_ << endl;
@@ -127 +144 @@
 void BRMeanSpecCalculator::DefineDataTable()
 {
+  cout << "(JEC) BRMeanSpecCalculator::DefineDataTable START" << endl;
   string dtfile="!"+outpath_+"/dtspec.fits";
   ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
@@ -136 +154 @@
     dtp_->AddFloatColumn(cnom);    
   }
+  for(int i=0; i<nchan; i++) {
+    for(int j=0;j<numberOfBands_;j++){
+      sprintf(cnom,"varnormb%d%d",i,j);
+      dtp_->AddFloatColumn(cnom);    
+    }
+  }
   /*
   for(int i=0; i<nchan; i++) {
@@ -143 +167 @@
   */
   //  xnt_=new double[nchan*2];  CHECK : faut-il reallouer ?
+  cout << "(JEC) creation tuple de " << nchan*(1+numberOfBands_) << " doubles " << endl;
+  xnt_=new double[nchan*(1+numberOfBands_)];
+  cout << "(JEC) BRMeanSpecCalculator::DefineDataTable END" << endl;
+}
+  
+/* --Methode-- */
+void BRMeanSpecCalculator::SetNumberOfBands(int numberOfBands, int ibandfirst, int ibandlast)
+{
+  if (numberOfBands < 1) numberOfBands = 1;
+  numberOfBands_ = numberOfBands;
+  if (ibandfirst < 0 )ibandfirst = 0;
+  if (ibandlast >= numberOfBands_ ) ibandlast = numberOfBands_-1;
+  ibandfirst_=ibandfirst; ibandlast_=ibandlast;
+
+  cout << "(JEC): SetNumberOfBands (END) : "
+       << numberOfBands_ << " "
+       << ibandfirst_ << " "
+       << ibandlast_ << endl;
+
 }
 
@@ -154 +197 @@
   cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << " ApplyGains=" 
        << ((fg_apply_gains_)?"true":"false") << endl;
-}
-
-static inline r_4 Zmod2(complex<r_4> z) 
-{ return (z.real()*z.real()+z.imag()*z.imag()); }
+  if( (spec_window_.SizeX()!= sgain_.NCols()) || (spec_window_.SizeY()!= sgain_.NRows()) ){
+    cout << " ReadGainFitsFile: BAD Gain Matrix sizes " << endl;
+    sgain_.Show();
+    spec_window_.Show();
+    throw ParmError("ReadGainFitsFile: BAD Gain Matrix sizes");
+  }
+}
+
+//JEC
+//static inline r_4 Zmod2(complex<r_4> z) 
+//{ return (z.real()*z.real()+z.imag()*z.imag()); }
 
 
@@ -169 +219 @@
   // des spectres le cas echeant ProcSpecWin()
 
+ 
   int_8 nbpaqdec = (int_8)totnbpaq_-(int_8)GetSpecWinOverlapSize();
   if ((nbpaqdec>0)&&(nbpaqdec%GetSpectraWindowSize()==0)) {
@@ -208 +259 @@
   if (fg_apply_gains_) {   // Application des gains, si demande 
     sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
-    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) 
-      (spec_window_(Range::all(), Range(i), Range(kz))).Div(sgain_.Row(i));
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
+      (spec_window_(Range::all(), Range(i), Range(kz)).CompactAllDimensions()).Div(sgain_.Row(i).CompactAllDimensions());
+    }
   }
 
@@ -231 +283 @@
     }
   }
-  nbtot_specwin_++;
-  return;
   // On appelle la routine de nettoyage qui doit flagger les mauvais paquets 
   FlagBadPackets(numpaqstart, numpaqend);
@@ -249 +299 @@
       // Calcul de spectres moyennes et variance 
       for(sa_size_t f=1; f<spec.Size(); f++)  {   // boucle sur les frequences 
-        spec(f) += spec_window_(f,i,kz);
+        spec(f)  += spec_window_(f,i,kz);
         sspec(f) += spec_window_(f,i,kz)*spec_window_(f,i,kz);  
       }
@@ -270 +320 @@
     for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  {
       double mean, sigma;
-      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
+      ////////BUG      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
       double variance=0.;
       variance=spec_window_(Range(1,Range::lastIndex()), Range(i), Range(kz)).Sum();
-      //      xnt_[i]=variance;
+      xnt_[i]=variance;
+      //Compute nomalized variance in bands freq. 
+      sa_size_t fMin;
+      sa_size_t fMax;
+      int bandW = spec_window_.SizeX()/numberOfBands_;
+      vector<double> varNomBinned(numberOfBands_);
+      for (sa_size_t j=0; j<numberOfBands_; j++){
+        fMin = j*bandW;
+        fMax =fMin+bandW-1;
+        varNomBinned[j]=spec_window_(Range(fMin,fMax), Range(i), Range(kz)).Sum();
+        varNomBinned[j]/=(r_4)bandW;
+        xnt_[spec_window_.SizeY()+i*numberOfBands_+j] =  varNomBinned[j];
+      }//eof
+
       clnflg_(i,kz)=0;
-      if (variance<varmin_) { clnflg_(i,kz)=1;  nbadpaq_[i]++; }
-      else if (variance>varmax_) { clnflg_(i,kz)=2;  nbadpaq_[i]++; }
+      for (sa_size_t j=ibandfirst_; j<=ibandlast_; j++){
+        //        cout << "(jec) var["<<j<<"] =" << varNomBinned[j]
+        //     << " min " << varmin_
+        //     << " max " << varmax_ << endl;
+        if(varNomBinned[j]<varmin_)      { clnflg_(i,kz)=10+j;  nbadpaq_[i]++; break;}
+        else if(varNomBinned[j]>varmax_) { clnflg_(i,kz)=100+j; nbadpaq_[i]++; break;}
+      }
+      //cout << "clnflg_("<<i<<","<<kz<<"): " << clnflg_(i,kz)  << endl;
     }
     if (dtp_)  dtp_->AddRow(xnt_);
@@ -343 +412 @@
   : BRBaseProcessor(memgr), fgsinglechannel_(fgsinglechannel), totnbfftpaq_(0)
 {
-  BRPaquet paq(memgr_.PaqSize());
+  uint_4 nb_octets_entrop = 0; //this value is valid for Dec. 2010 data at Nancay
+  const char* venvp = NULL; 
+  venvp=getenv("BRANA_NBYTECUT");
+  if (venvp!=NULL){
+    nb_octets_entrop = atoi(venvp);
+    cout << "BRFFTCalculator  : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
+  }
+
+  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
+  //JEC END
+  //  BRPaquet paq(memgr_.PaqSize());
   setNameId("FFTCalc",2);
   ffts_.SetInDataSize((fgsinglechannel_)?paq.DataSize():paq.DataSize()/2);
@@ -633 +712 @@
         }
         for(sa_size_t j=0; j<spectreV1.Size(); j++) 
-          spectreV1(j) += Zmod2(cfour1(j));
+          spectreV1(j) += BRMeanSpecCalculator::Zmod2(cfour1(j));
         memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
         if (fgtimfreq) {   // Remplissage tableau temps-frequence
           for(sa_size_t c=1; c<timfreqV1.NCols(); c++) {
             for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) 
-              timfreqV1(nzm, c) += Zmod2(cfour1(j));
+              timfreqV1(nzm, c) += BRMeanSpecCalculator::Zmod2(cfour1(j));
           }
         }
         for(sa_size_t j=0; j<spectreV2.Size(); j++) 
-          spectreV2(j) += Zmod2(cfour2(j));  // Zmod2(zp2[j]); 
+          spectreV2(j) += BRMeanSpecCalculator::Zmod2(cfour2(j));  // BRMeanSpecCalculator::Zmod2(zp2[j]); 
         memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
         if (fgtimfreq) {   // Remplissage tableau temps-frequence
           for(sa_size_t c=1; c<timfreqV2.NCols(); c++) {
             for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) 
-              timfreqV2(nzm,c) += Zmod2(cfour2(j));
+              timfreqV2(nzm,c) += BRMeanSpecCalculator::Zmod2(cfour2(j));
           }
         }
@@ -667 +746 @@
         if (nmnt==0)  { xnt[0]=paq.FrameCounter();  xnt[1]=paq.TimeTag(); }
         for(sa_size_t j=2700; j<2800; j++) { 
-          ms1 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j)); 
+          ms1 += BRMeanSpecCalculator::Zmod2(cfour1(j)); ms2 += BRMeanSpecCalculator::Zmod2(cfour2(j)); 
           complex<r_4> zvis =  cfour1(j)*conj(cfour2(j));
-          ms12 += Zmod2(zvis);   ms12re += zvis.real();  ms12im += zvis.imag();
+          ms12 += BRMeanSpecCalculator::Zmod2(zvis);   ms12re += zvis.real();  ms12im += zvis.imag();
           ms12phi+= atan2(zvis.imag(),zvis.real());
         } 
@@ -742 +821 @@
           int_4 sspvmaxidx[3] = {-1,-1,-1};
           for(int jji=1;jji<visiV12.Size()-1;jji++) {
-            r_4 zmv2 = Zmod2(visiV12(jji)); 
+            r_4 zmv2 = BRMeanSpecCalculator::Zmod2(visiV12(jji)); 
             if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; }
           }