Index: /trunk/AddOn/TAcq/branap.cc
===================================================================
--- /trunk/AddOn/TAcq/branap.cc	(revision 3942)
+++ /trunk/AddOn/TAcq/branap.cc	(revision 3943)
@@ -26,5 +26,5 @@
   nzones_=nzon;
   npaqinzone_=npaqz;
-  fgdatafft_=false;  fgsinglechannel_=false;
+  fgdatafft_=false;  fgsinglechannel_=false; fgforcesinglechan_=false;
   prtlevel_=0;
   prtmodulo_=50000;
@@ -37,4 +37,9 @@
 
   vmin_=0.;  vmax_=9e99;
+  nbands_=0; bandfirst_ = bandlast_ = 0;
+  fgdt_ = false;
+
+  fgfreqfilter_=false; //JEC 1/2/11
+
   gainfile_="";
 }
@@ -94,8 +99,24 @@
       ka++;
     }
+    else if (strcmp(arg[ka],"-singlechan")==0) {
+      fgsinglechannel_=true; fgforcesinglechan_ = true; 
+      ka++;
+    }
     else if (strcmp(arg[ka],"-varcut")==0) {
       sscanf(arg[ka+1],"%lg,%lg",&vmin_,&vmax_);
       ka+=2;
     }
+    else if (strcmp(arg[ka],"-nband")==0) {
+      sscanf(arg[ka+1],"%d,%d,%d",&nbands_,&bandfirst_,&bandlast_);
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-filldt")==0) {
+      fgdt_=true;
+      ka++;
+    }
+    else if (strcmp(arg[ka],"-freqfilter")==0) {
+      fgfreqfilter_=true;
+      ka++;
+    }    
     else if (strcmp(arg[ka],"-gain")==0) {
       gainfile_=arg[ka+1];
@@ -147,6 +168,8 @@
   cout << " Usage:  prgname [-act ACT] [-out OutPath] [-nmean NMean] [-zones NZones,nPaqinZone] \n" 
        << "                 [-nbloc NBloc] [-freq NumFreqMin,NumFreqMax,NBinFreq] \n" 
-       << "                 [-prt lev,modulo] [-nvcal n] [-nthr n] [-nosfc] \n" 
-       << "                 [-gain filename] [-varcut min,max] [-tspwin wsz,extsz,nfiles] \n" 
+       << "                 [-prt lev,modulo] [-nvcal n] [-nthr n] [-nosfc] [-singlechan]\n"
+       << "                 [-freqfilter]\n"
+       << "                 [-gain filename] [-varcut min,max] [-nband nband,first,last] \n"
+       << "                 [-filldt] [-tspwin wsz,extsz,nfiles] \n" 
        << "         -in Imin,Imax,Istep InPath FiberList [InPath2 FiberList2 InPath3 FiberList3 ...] \n" << endl; 
   if (fgshort) {
@@ -168,5 +191,9 @@
        << " -nthr n : number of threads for parallel execution in BRVisibilityCalculator (default=1) \n" 
        << " -nosfc : Don't force reading with SAME FrameCounter \n" 
+       << " -singlechan : Force one channel per fiber \n"
        << " -varcut min,max : min-max cut on variance \n" 
+       << " -nband nband,first,last: numbers of freq. bands and first and last bands used for cuts \n"
+       << " -filldt : force data table filling \n"
+       << " -freqfilter : force median filtering on the frequencies \n"
        << " -gain filename : spectral response fits file name \n" 
        << " -tspwin wsz,extsz,nfiles : spectra time (paquet no) window (ex: -tspwin 120,4,5) \n"
@@ -185,5 +212,7 @@
   char flnm[1024];
   sprintf(flnm,"%s/signal%d.fits",dirlist_[0].c_str(),imin_);
-  return DecodeMiniFitsHeader(flnm,paqsize_, npaqinfile_,fgdatafft_, fgsinglechannel_); 
+  int rc = DecodeMiniFitsHeader(flnm,paqsize_, npaqinfile_,fgdatafft_, fgsinglechannel_); 
+  if(fgforcesinglechan_) fgsinglechannel_= true;
+  return rc;
 }
 
@@ -199,5 +228,8 @@
   cout << " Action=" << action_ << "  NMean=" << nmean_ << " NBloc=" << nbloc_ << endl;
   cout << " FreqMin= " << freqmin_ << " FreqMax= " << freqmax_ << " NBinFreq= " << nbinfreq_ << endl;
-  cout << " GainFileName=" << gainfile_ << " VarianceMin= " << vmin_ << " VarianceMax= " << vmax_ << endl;
+  cout<<  ((fgdt_)?" Fill DadaTable ":" NO DataTable") << endl;
+  cout << " GainFileName=" << gainfile_ << " Variance: Min= " << vmin_ << " Max= " << vmax_ 
+       << " Bands: N=" <<  nbands_ << " First=" << bandfirst_ << " Last=" << bandlast_
+       << endl;
   cout << " Spectra TimeWindow (Nb.Paquets) Size=" << spec_win_sz_ << " ExtensionSize=" << spw_ext_sz_ 
        << " MaxNbFile=" << nbmax_specwfiles_ << endl; 
Index: /trunk/AddOn/TAcq/branap.h
===================================================================
--- /trunk/AddOn/TAcq/branap.h	(revision 3942)
+++ /trunk/AddOn/TAcq/branap.h	(revision 3943)
@@ -41,4 +41,5 @@
   bool fgdatafft_;  // true -> donnee FFT (provenant du firmware FFT)
   bool fgsinglechannel_; // true -> un seul canal par fibre (par defaut=2 canaux/fibres)
+  bool fgforcesinglechan_; // true force par la ligne de commande
   int prtlevel_;
   long int prtmodulo_; 
@@ -51,4 +52,9 @@
 
   double vmin_,vmax_;   // coupure min-max sur la variance 
+  int nbands_,bandfirst_,bandlast_; //bandes en frequences pour mean, variance et cut
+  bool fgdt_; //if true fill datatable
+
+  bool fgfreqfilter_; // if true force median filtering on the frequencies
+
   string gainfile_;     // nom du fichier de gain (reponse spectrale)
 };
Index: /trunk/AddOn/TAcq/brproc.cc
===================================================================
--- /trunk/AddOn/TAcq/brproc.cc	(revision 3942)
+++ /trunk/AddOn/TAcq/brproc.cc	(revision 3943)
@@ -43,8 +43,18 @@
     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,13 +80,20 @@
   }
 
+  //
+  
+
   // 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,8 +101,8 @@
 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,4 +144,5 @@
 void BRMeanSpecCalculator::DefineDataTable()
 {
+  cout << "(JEC) BRMeanSpecCalculator::DefineDataTable START" << endl;
   string dtfile="!"+outpath_+"/dtspec.fits";
   ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
@@ -136,4 +154,10 @@
     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,4 +167,23 @@
   */
   //  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,8 +197,15 @@
   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,4 +219,5 @@
   // des spectres le cas echeant ProcSpecWin()
 
+ 
   int_8 nbpaqdec = (int_8)totnbpaq_-(int_8)GetSpecWinOverlapSize();
   if ((nbpaqdec>0)&&(nbpaqdec%GetSpectraWindowSize()==0)) {
@@ -208,6 +259,7 @@
   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,6 +283,4 @@
     }
   }
-  nbtot_specwin_++;
-  return;
   // On appelle la routine de nettoyage qui doit flagger les mauvais paquets 
   FlagBadPackets(numpaqstart, numpaqend);
@@ -249,5 +299,5 @@
       // 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,11 +320,30 @@
     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,5 +412,15 @@
   : 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,19 +712,19 @@
 	}
 	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,7 +746,7 @@
         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,5 +821,5 @@
 	  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; }
 	  }
Index: /trunk/AddOn/TAcq/brproc.h
===================================================================
--- /trunk/AddOn/TAcq/brproc.h	(revision 3942)
+++ /trunk/AddOn/TAcq/brproc.h	(revision 3943)
@@ -37,5 +37,5 @@
 //---------------------------------------------------------------------
 
-class BRMeanSpecCalculator : public BRBaseProcessor
+class BRMeanSpecCalculator : public BRBaseProcessor 
 {
 public:
@@ -43,4 +43,7 @@
 		       bool fgdatafft=true, bool fgsinglechan=false);
   virtual ~BRMeanSpecCalculator(); // Fait le nettoyage final , sauvegarde des objets ...
+
+ //JEC 19/1/11 
+  static inline r_4 Zmod2(complex<r_4> z){ return (z.real()*z.real()+z.imag()*z.imag()); }
 
 // Definit la taille de la fenetre temporelle (winsz = nb de paquets) des spectres a traiter
@@ -55,7 +58,8 @@
 // Retourne l'index (selon Z) dans le tableau des spectres en fonction du numero de paquet
   inline sa_size_t PaqNumToArrayIndex(uint_8 numpaq)  { return numpaq%spec_window_.SizeZ() ; }
+  inline int GetNumberOfBands() { return numberOfBands_; }
 
 // Pour definir le nombre maximum de fichiers de fenetres de spectres a sauver sur fichiers
-  inline void SetMaxNbSepcWinFiles(uint_4 nmax=0)  {  nmaxfiles_specw_=nmax; }
+  inline void SetMaxNbSpecWinFiles(uint_4 nmax=0)  {  nmaxfiles_specw_=nmax; }
 
 // Pour definir le datatable a remplir - pas de DataTable rempli si pas appele 
@@ -64,4 +68,9 @@
   inline void SetVarianceLimits(double vmin=0., double vmax=9.e99)
     { varmin_=vmin;  varmax_=vmax; }
+  
+  //POur definir le nbre de bandes de freequences pour calcul mean, variance, cut
+  //ibandlast incluse
+  void SetNumberOfBands(int numberOfBands = 1, int ibandfirst = 0, int ibandlast = 99999);
+
 
 // lecture des gains en fonction de la frequence depuis un fichier FITS
@@ -77,5 +86,5 @@
   virtual void SaveSpectraWindow();  // Pour sauver les spectres de la fenetre temporel ds un fichier 
 
-
+   
   uint_4 nmean_;  // Nombre de spectres pour le calcul des moyennes 
   string outpath_;  // directory pour fichiers de sortie
@@ -97,4 +106,5 @@
   uint_4 nbtot_specwin_;    // nombre total de fenetre de spectres traites 
   uint_4 nmaxfiles_specw_;   // nombre maximum de fichier de fenetre_spectre 
+
   
   TMatrix< r_4 > sgain_;   // gain en fonction de la frequence  
@@ -102,4 +112,7 @@
 
   double varmin_, varmax_;  // limites pour la coupure sur la variance
+  int ibandfirst_, ibandlast_; //limites pour les numero de bandes a utiliser pour les cut en varaince  
+  int numberOfBands_; // nombre de bandes de frequences pour calcul mean, variance et cut
+
   vector<uint_8> nbadpaq_;
 
Index: /trunk/AddOn/TAcq/brprocGain.cc
===================================================================
--- /trunk/AddOn/TAcq/brprocGain.cc	(revision 3943)
+++ /trunk/AddOn/TAcq/brprocGain.cc	(revision 3943)
@@ -0,0 +1,196 @@
+//----------------------------------------------------------------
+// Projet BAORadio - (C) LAL/IRFU  2008-2010
+// Classes de threads de traitement pour BAORadio 
+//----------------------------------------------------------------
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <fstream>
+#include <signal.h>
+
+#include "pexceptions.h"
+#include "tvector.h"
+#include "ntuple.h"
+#include "datatable.h"
+#include "histos.h"
+#include "fioarr.h"
+#include "matharr.h"
+#include "timestamp.h"
+#include "ctimer.h"
+#include "fftpserver.h"
+#include "fitsarrhand.h"
+
+#include "FFTW/fftw3.h"
+
+
+#include "pciewrap.h"
+#include "brpaqu.h"
+
+#include "brprocGain.h"
+
+
+
+
+//---------------------------------------------------------------------
+// Classe de traitement de spectres - 
+// Calcul de spectres de gain
+//---------------------------------------------------------------------
+/* --Methode-- */
+BRGainCalculator::BRGainCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, 
+					   bool fgdatafft, bool fgsinglechan)
+  : BRMeanSpecCalculator(memgr,outpath,nmean,fgdatafft,fgsinglechan), forceMedianFreqFilter_(false),nbwin4mean_(nmean)
+{
+
+  //  cout << "(JEC): BRGainCalculator Ctor " << endl;
+
+  
+  setNameId("gainCalc",1);
+
+  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 << "BRGainCalculator : BRANA_NBYTECUT : " << nb_octets_entrop << endl;
+  }
+
+  BRPaquet paq(memgr_.PaqSize()-nb_octets_entrop);
+
+  if (fgsinglechannel_) {
+    medfiltspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2);     
+  }
+  else {
+    medfiltspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4);     
+  }
+  //init
+  medfiltspecmtx_ = (r_4)(0.);
+  nummedianfile_ = 0;
+  nbtot_specwin_ = 0;
+
+}
+
+
+/* --Methode-- */
+BRGainCalculator::~BRGainCalculator()
+{
+  cout << "(JEC): BRGainCalculator Dtor:  nbtot_specwin = " << nbtot_specwin_
+       << "nbwin4mean = " << nbwin4mean_ << endl;
+  if(nbtot_specwin_>nbwin4mean_*0.5)SaveMedianSpectra();
+}
+
+
+
+/* --Methode-- */
+void BRGainCalculator::ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend)
+{
+
+  cout << "(JEC): BRGainCalculator ProcSpecWin: "<<   nbtot_specwin_ 
+       << " nbwin4mean = " << nbwin4mean_
+       << endl;
+  
+  //  cout<< "(Spec_window info) Start" << endl;
+  //spec_window_.Show();
+  //cout<< "(Spec_window info) End" << endl;
+
+  //JEC 13/12/10 save the filtered spectra
+  if ( (nbtot_specwin_>0)&&(nbtot_specwin_%nbwin4mean_==0) ) SaveMedianSpectra();
+
+
+  //DBG  cout << "BRGainCalculator::ProcSpecWin()/Debug: numpaqstart=" << numpaqstart 
+  //DBG     << " numpaqend=" << numpaqend << endl;
+
+
+  //Implement the median filter at a fixed freq. but on all paquets of the window
+  //then apply a median filter on the frequencies
+  for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  { //loop on channel (voie)
+    TVector<r_4> medianFilterPaq(spec_window_.SizeX());
+    //we keep freq. = 0 for compatibility
+    for (sa_size_t f=0; f<spec_window_.SizeX(); f++) {
+      std::vector<r_4> val1;
+      TVector<r_4> tval1(spec_window_(Range(f),Range(i),Range::all()).CompactAllDimensions());
+      tval1.FillTo(val1);
+      std::nth_element(val1.begin(),
+		       val1.begin()+val1.size()/2,
+		       val1.end());
+      medianFilterPaq(f) = *(val1.begin()+val1.size()/2);
+    }//eo loop on frequencies
+
+
+
+    if(forceMedianFreqFilter_) {
+      //perform a median filter on the frequencies
+      TVector<r_4> medianFilterFreq = medfiltspecmtx_.Row(i);
+      sa_size_t hww = 50; //half freq window for filtering
+      sa_size_t fMin = 0.;
+      sa_size_t fMax = spec_window_.SizeX()-1;
+      for (sa_size_t f=0; f<spec_window_.SizeX(); f++) {
+	sa_size_t first = (f-hww >= fMin) ? f-hww : fMin;
+	sa_size_t last = (f+hww <= fMax) ? f+hww : fMax;
+	std::vector<r_4> val2;
+	medianFilterPaq(Range(first,last)).FillTo(val2);
+	std::nth_element(val2.begin(),
+			 val2.begin()+val2.size()/2,
+			 val2.end());
+	medianFilterFreq(f) += *(val2.begin()+val2.size()/2); //Notice the "+" for later take the mean
+      }//eo loop on frequencies
+    } else {
+      //      cout<< "(JEC) ProcSpecWin debug medfiltspecmtx/medianFilterPaq "<< endl;
+      //      medfiltspecmtx_.Row(i).Show();
+      //      medianFilterPaq.Transpose().Show();
+      medfiltspecmtx_.Row(i) = medianFilterPaq.Transpose();
+    }
+  }//eo loop on channels
+
+  //We always end here with a medfiltspecmtx Matrix Row=[0,NCha-1] x Col=[0,fMax]
+
+  
+  if (nbtot_specwin_<nmaxfiles_specw_)  SaveSpectraWindow();
+  
+
+  nbtot_specwin_++;
+  return;
+}
+
+
+/* --Methode-- */
+//JEC 13/12/10 Start
+void BRGainCalculator::SaveMedianSpectra()
+{
+
+  cout << "(JEC): BRGainCalculator   SaveMedianSpectra ("<<nummedianfile_<<") nbwin = "<< nbwin4mean_  
+       << " nbtot_specwin = " <<  nbtot_specwin_
+       << endl;
+  
+    for(sa_size_t ir=0; ir< medfiltspecmtx_.NRows(); ir++){
+      char buff[32];
+      sprintf(buff,"NPAQSUM_%d",(int)ir);
+      medfiltspecmtx_.Info()["NPAQSUM"] = nbtot_specwin_;
+      medfiltspecmtx_.Info()[buff] = nbtot_specwin_;
+      if ( nbtot_specwin_ > 0) {
+	medfiltspecmtx_.Row(ir)  /= (r_4)nbtot_specwin_;
+      }
+    }//eo loop on channels
+    char nfile[64];
+    string flnm;
+    {
+      sprintf(nfile,"medfiltmtx%d.fits",nummedianfile_);
+      flnm="!"+outpath_+nfile;
+      FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
+      fos << medfiltspecmtx_;
+    }
+    cout << nummedianfile_ << "-BRGainCalculator::ProcSpecWin() save filtered spectra matrix in " 
+	 << flnm << endl;
+
+    //increment the file index
+    nummedianfile_++;
+    
+    //reset the matirx
+    medfiltspecmtx_ = (r_4)(0.);
+
+    //reset counter
+    nbtot_specwin_ = 0;
+
+    return;
+}
+//JEC 13/12/10 End
Index: /trunk/AddOn/TAcq/brprocGain.h
===================================================================
--- /trunk/AddOn/TAcq/brprocGain.h	(revision 3943)
+++ /trunk/AddOn/TAcq/brprocGain.h	(revision 3943)
@@ -0,0 +1,47 @@
+#ifndef  BRPROCGAIN_H_SEEN
+#define  BRPROCGAIN_H_SEEN
+
+//----------------------------------------------------------------
+// Projet BAORadio - (C) LAL/IRFU  2008-2011
+// Classes de threads de traitememt donnees BAORadio 
+//  Principales classes utilisables 
+//  - BRGainCalculator (donnees FFT ou RAW apres BRFFTCalculator)
+//  - BRFFTCalculator (donnees RAW) 
+//----------------------------------------------------------------
+
+#include "brproc.h"
+
+using namespace std;
+
+ 
+//---------------------------------------------------------------------
+// Classe de traitement simple pour obtenir des spectres filtres des RFIs 
+// que l'on appelle spectre de gain / voie
+// a partir de donnees de type FFT (firmware FFT)
+//---------------------------------------------------------------------
+class BRGainCalculator : public BRMeanSpecCalculator
+{
+public:
+  BRGainCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean=1000, 
+  		       bool fgdatafft=true, bool fgsinglechan=false);
+  virtual ~BRGainCalculator(); // Fait le nettoyage final , sauvegarde des objets ...
+
+  void SetOnMedianFreqFilt(){forceMedianFreqFilter_ = true;}  
+
+
+protected:
+
+// appele par Process() pour traiter un groupe de paquets correspondant a une fenetre en temps
+  virtual void ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend);  // appele fenetre par fenetre  
+
+  virtual void SaveMedianSpectra();  // Pour sauver les spectres median dans un fichier
+
+
+  bool forceMedianFreqFilter_; // if true force Median frequency filtering
+
+  TMatrix< r_4 > medfiltspecmtx_;    // Matrice des spectres filtres 
+  uint_8 nbwin4mean_;    // number of windows to perform the mean value of the filtered spectra
+  uint_4 nummedianfile_; //index of the filtered spectra file 
+
+}; 
+#endif
Index: /trunk/AddOn/TAcq/makefile
===================================================================
--- /trunk/AddOn/TAcq/makefile	(revision 3942)
+++ /trunk/AddOn/TAcq/makefile	(revision 3943)
@@ -5,5 +5,5 @@
 
 PKGHLIST = bracqvers.h branap.h brparam.h racqumem.h brpaqu.h brtypes.h \
-           brdiskw.h brfitsrd.h minifits.h brviscalc.h brbaseproc.h brproc.h \
+           brdiskw.h brfitsrd.h minifits.h brviscalc.h brbaseproc.h brproc.h brprocGain.h\
            mfacq.h  tmtacq.h racqueth.h  racquproc.h racqurw.h pciewrap.h  \
            swrapsock.h burawstream.h sockrawstream.h tmtacq.h
@@ -11,5 +11,5 @@
 PKGOLIST = $(OBJ)branap.o  $(OBJ)brparam.o $(OBJ)racqumem.o $(OBJ)brpaqu.o \
            $(OBJ)brfitsrd.o $(OBJ)brdiskw.o  $(OBJ)minifits.o \
-           $(OBJ)brbaseproc.o $(OBJ)brproc.o $(OBJ)brviscalc.o \
+           $(OBJ)brbaseproc.o $(OBJ)brproc.o $(OBJ)brprocGain.o $(OBJ)brviscalc.o \
            $(OBJ)racqueth.o $(OBJ)racquproc.o $(OBJ)racqurw.o  $(OBJ)pciewrap.o \
            $(OBJ)swrapsock.o $(OBJ)sockrawstream.o $(OBJ)burawstream.o
@@ -115,4 +115,7 @@
 	$(CXXCOMPILE) -c  -o $(OBJ)brviscalc.o brviscalc.cc 
 
+$(OBJ)brprocGain.o : brprocGain.cc $(PKGHLIST)
+	$(CXXCOMPILE) -c  -o $(OBJ)brprocGain.o brprocGain.cc 
+
 $(OBJ)brproc.o : brproc.cc $(PKGHLIST)
 	$(CXXCOMPILE) -c  -o $(OBJ)brproc.o brproc.cc 
Index: /trunk/AddOn/TAcq/specmfib.cc
===================================================================
--- /trunk/AddOn/TAcq/specmfib.cc	(revision 3942)
+++ /trunk/AddOn/TAcq/specmfib.cc	(revision 3943)
@@ -26,4 +26,9 @@
 #include "brfitsrd.h"
 #include "brproc.h"
+
+//JEC 19/1/11 START
+#include "brprocGain.h"
+//JEC END
+
 #include "brdiskw.h"
 
@@ -50,5 +55,9 @@
     if (rc) return rc;
     par.Print(cout);
-    if ((par.action_!="cube3d")&&(par.action_!="mspec")&&(par.action_!="bproc")) {
+    if ((par.action_!="cube3d")
+	&&(par.action_!="mspec")
+	&&(par.action_!="gain") //JEC 19/1/11 add gain action
+	&&(par.action_!="bproc") 
+	) {
       cout << " !!! specmfib.cc  BAD action = " << par.action_ << " possible values: mspec,cube3d,bproc" << endl;
       return 5;
@@ -73,14 +82,34 @@
     BRMeanSpecCalculator procms(mmgr, par.outpath_, par.nmean_, par.fgdatafft_, par.fgsinglechannel_);
     procms.SetVarianceLimits(par.vmin_, par.vmax_);
-    if (par.gainfile_.length()>0) procms.ReadGainFitsFile(par.gainfile_);
+    if (par.nbands_>0) procms.SetNumberOfBands(par.nbands_,par.bandfirst_,par.bandlast_);
+    //JEC 27/1/11 see below    if (par.gainfile_.length()>0) procms.ReadGainFitsFile(par.gainfile_);
     procms.SetPrintLevel(par.prtlevel_,par.prtmodulo_);
-    // Il faut appeler procms.DefineDataTable() si on veut remplir le DataTable
+    if (par.fgdt_) procms.DefineDataTable();
+
+    BRGainCalculator* procgain_p=NULL;
+    if (par.action_=="gain") {
+      procgain_p = 
+	new BRGainCalculator(mmgr, par.outpath_, par.nmean_, par.fgdatafft_, par.fgsinglechannel_);
+      procgain_p->SetPrintLevel(par.prtlevel_,par.prtmodulo_);
+    }
+
 
     BRFFTCalculator procfft(mmgr, par.fgsinglechannel_);
-    if (!par.fgdatafft_) procms.SetMemZAction(MemZA_ProcB);
+    if (!par.fgdatafft_) { 
+      procms.SetMemZAction(MemZA_ProcB);
+      if(procgain_p) procgain_p->SetMemZAction(MemZA_ProcB);
+    }
+
     if (par.spec_win_sz_>0) { 
       procms.SetSpectraWindowSize(par.spec_win_sz_, par.spw_ext_sz_);
-      procms.SetMaxNbSepcWinFiles(par.nbmax_specwfiles_);
+      procms.SetMaxNbSpecWinFiles(par.nbmax_specwfiles_);
+      if(procgain_p){
+	procgain_p->SetSpectraWindowSize(par.spec_win_sz_, par.spw_ext_sz_);
+	procgain_p->SetMaxNbSpecWinFiles(par.nbmax_specwfiles_);
+      }
     }
+    //JEC 27/1/11 should be done after SetSpectraWindowSize to perform x-checking
+    if (par.gainfile_.length()>0) procms.ReadGainFitsFile(par.gainfile_);
+
     FitsCubeWriter wrt(mmgr, par.outpath_, par.nbloc_);
 
@@ -92,11 +121,14 @@
     }
 
-    cout << " specmfib: Starting threads (reader meanSpecCalculator ...) ... " << endl;
+    cout << " specmfib: Starting threads (reader + "<< par.action_ << ")" << endl;
     reader.start();
     if (par.action_ == "cube3d")   wrt.start();
     else {  // Calcul spectre moyenne 
       if (!par.fgdatafft_)  procfft.start(); 
+      //JEC 19/1/11 START add the gain process alternative
       if (par.action_=="bproc")  bproc_p->start();
+      else if (par.action_=="gain") procgain_p->start();
       else procms.start();
+      //JEC END
     }
     usleep(200000);
@@ -105,9 +137,14 @@
     else {
       if (!par.fgdatafft_)  procfft.join(); 
+
+      //JEC 19/1/11 START add the gain process alternative
       if (par.action_=="bproc")  bproc_p->join();
+      else if (par.action_=="gain") procgain_p->join();
       else procms.join();
+      //JEC END
     }
     mmgr.Print(cout);
     if (bproc_p) delete bproc_p;
+    if (procgain_p) delete procgain_p;
     cout << resu ; 
   }
