Index: /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc
===================================================================
--- /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc	(revision 558)
+++ /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc	(revision 558)
@@ -0,0 +1,546 @@
+// Utilisation de SOPHYA pour faciliter les tests ...
+#include <regex.h>
+//#include <regexp.h>
+#include <stdio.h>
+
+#include "sopnamsp.h"
+#include "machdefs.h"
+
+#include <stdlib.h>
+#include <dirent.h>
+#include <matharr.h>
+
+// include standard c/c++
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <map>
+#include <functional>
+#include <algorithm>
+#include <numeric>
+#include <list>
+#include <exception>
+
+// include sophya mesure ressource CPU/memoire ...
+#include "resusage.h"
+#include "ctimer.h"
+#include "timing.h"
+#include "timestamp.h"
+#include "strutilxx.h"
+#include "ntuple.h"
+#include "fioarr.h"
+#include "tarrinit.h"
+#include "histinit.h"
+#include "fitsioserver.h"
+#include "fiosinit.h"
+#include "ppersist.h"
+
+//-----------------------------------------------
+const sa_size_t NUMBER_OF_CHANNELS = 2;
+const sa_size_t NUMBER_OF_FREQ     = 8192;
+const sa_size_t TOTAL_NUM_CYCLES   = 500;
+const r_4    LOWER_FREQUENCY = 1250.0; //MHz
+const r_4    TOTAL_BANDWIDTH = 250.0; //MHz
+//-----------------------------------------------
+//Input parameters
+struct Param {
+  int debuglev_;      //debug
+  string inPath_;     //root directory of the input files
+  string outPath_;    //output files are located here
+  string sourceName_; //source name & subdirectory of the input files
+  string ppfFile_;    //generic name of the input files
+  int nSliceInFreq_;  //used by reduceSpectra() fnc
+} para;
+//--------------------------------------------------------------
+//Utility functions
+
+sa_size_t freqTochan(r_4 f){
+  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
+}
+//---------
+  class median_of_empty_list_exception:public std::exception{
+    virtual const char* what() const throw() {
+    return "Attempt to take the median of an empty list of numbers.  "
+      "The median of an empty list is undefined.";
+  }
+  };
+  template<class RandAccessIter>
+    double median(RandAccessIter begin, RandAccessIter end) 
+    throw(median_of_empty_list_exception){
+    if(begin == end){ throw median_of_empty_list_exception(); }
+    std::size_t size = end - begin;
+    std::size_t middleIdx = size/2;
+    RandAccessIter target = begin + middleIdx;
+    std::nth_element(begin, target, end);
+    
+    if(size % 2 != 0){ //Odd number of elements
+      return *target;
+    }else{            //Even number of elements
+      double a = *target;
+      RandAccessIter targetNeighbor= target-1;
+      std::nth_element(begin, targetNeighbor, end);
+      return (a+*targetNeighbor)/2.0;
+    }
+  }
+
+//--------------------------------------------------------------
+char *regexp (const char *string, const char *patrn, int *begin, int *end) {    
+        int i, w=0, len;                 
+        char *word = NULL;
+        regex_t rgT;
+        regmatch_t match;
+        regcomp(&rgT,patrn,REG_EXTENDED);
+        if ((regexec(&rgT,string,1,&match,0)) == 0) {
+                *begin = (int)match.rm_so;
+                *end = (int)match.rm_eo;
+                len = *end-*begin;
+                word=(char*)malloc(len+1);
+                for (i=*begin; i<*end; i++) {
+                        word[w] = string[i];
+                        w++; }
+                word[w]=0;
+        }
+        regfree(&rgT);
+        return word;
+}
+//------- 
+sa_size_t round_sa(r_4 r) {
+  return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
+}
+//-----
+string StringToLower(string strToConvert){
+  //change each element of the string to lower case
+  for(unsigned int i=0;i<strToConvert.length();i++) {
+    strToConvert[i] = tolower(strToConvert[i]);
+  }
+  return strToConvert;//return the converted string
+}
+//-----
+bool stringCompare( const string &left, const string &right ){
+   if( left.size() < right.size() )
+      return true;
+   for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
+      if( tolower( *lit ) < tolower( *rit ) )
+         return true;
+      else if( tolower( *lit ) > tolower( *rit ) )
+         return false;
+   return false;
+}
+//-----
+list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
+  list<string> theList;
+
+
+  DIR *dip;
+  struct dirent *dit;
+  string msg;  string fileName;
+  string fullFileName;
+  size_t found;
+
+  if ((dip=opendir(dir.c_str())) == NULL ) {
+    msg = "opendir failed on directory "+dir;
+    throw msg;
+  }
+  while ( (dit = readdir(dip)) != NULL ) {
+    fileName = dit->d_name;
+    found=fileName.find(filePettern);
+    if (found != string::npos) {
+      fullFileName = dir + "/";
+      fullFileName += fileName;
+      theList.push_back(fullFileName);
+    }
+  }//eo while
+  if (closedir(dip) == -1) {
+    msg = "closedir failed on directory "+dir;
+    throw msg;
+  }
+  
+  theList.sort(stringCompare);
+
+  return theList;
+
+}
+//
+class  StringMatch : public unary_function<string,bool> {
+public:
+  StringMatch(const string& pattern): pattern_(pattern){}
+  bool operator()(const string& aStr) const {
+
+
+    int b,e;
+    regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
+
+//     cout << "investigate " << aStr << " to find " << pattern_ 
+// 	 << "[" <<b<<","<<e<<"]" 
+// 	 << endl;
+
+    
+    if (b != 0) return false;
+    if (e != (int)aStr.size()) return false;
+    return true;
+
+  }
+private:
+  string pattern_;
+};
+//-------------------------------------------------------
+//Rebin in frequence + compute mean and sigma
+void reduceSpectra(const TMatrix<r_4>& specMtxInPut, 
+		    TMatrix<r_4>& meanMtx, 
+		    TMatrix<r_4>& sigmaMtx) {
+  sa_size_t nSliceFreq = para.nSliceInFreq_;
+  sa_size_t deltaFreq =  NUMBER_OF_FREQ/nSliceFreq;
+  for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
+    sa_size_t freqLow= iSlice*deltaFreq;
+    sa_size_t freqHigh= freqLow + deltaFreq -1;
+    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+      TVector<r_4> reducedRow;
+      reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
+      double mean; 
+      double sigma;
+      MeanSigma(reducedRow,mean,sigma);
+      meanMtx(iCh,iSlice) = mean;
+      sigmaMtx(iCh,iSlice) = sigma;
+    }//eo loop on channels
+  }//eo loop on slices
+}
+//-------------------------------------------------------
+//Compute the mean of Raw spectra and also the mean/sigma of rebinned spectra 
+//Used like:
+//
+void meanOnCycles() throw(string) {
+  list<string> listOfFiles;
+  string directoryName;
+  directoryName = para.inPath_ + "/" + para.sourceName_;
+
+  //Make the listing of the directory
+  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
+  
+  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
+  iFileEnd = listOfFiles.end();
+  
+  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
+  TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+  uint_4 nSpectra=0;
+  //Loop on files
+  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
+    if (para.debuglev_>90){
+      cout << "load file <" << *iFile << ">" << endl;
+    }
+    PInPersist fin(*iFile);
+    vector<string> vec = fin.GetNameTags();
+    list<string> listOfSpectra;
+    //Keep only required PPF objects
+    std::remove_copy_if(
+			vec.begin(), vec.end(), back_inserter(listOfSpectra),
+			not1(match)
+			);
+   
+    iSpecEnd = listOfSpectra.end();
+    listOfSpectra.sort(stringCompare);
+    //Loop of spectra matrix
+    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
+      if (para.debuglev_>90){
+	cout << " spactra <" << *iSpec << ">" << endl;
+      }
+      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      fin.GetObject(aSpec,*iSpec);
+      //How to see if the GetObject is ok?? Ask Reza
+      nSpectra++;
+      meanOfSpectra+=aSpec;
+    }//eo loop on spectra in a file
+  }//eo loop on files
+  
+  //Normalisation
+  if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
+
+  //Compute the reduced version of the mean and sigma
+  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
+
+  {//Save the result
+    stringstream tmp;
+    tmp << nSpectra;
+    string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+    cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
+    POutPersist fos(fileName);
+
+    string tag = "mean";
+    fos << PPFNameTag(tag) << meanOfSpectra;
+    tag = "meanred";
+    fos << PPFNameTag(tag) << meanRedMtx;
+    tag = "sigmared";
+    fos << PPFNameTag(tag) << sigmaRedMtx;
+  }
+}
+//-------------------------------------------------------
+//Compute the median of Raw spectra and also the mean/sigma of rebinned spectra 
+//Used like:
+//
+void medianOnCycles() throw(string) {
+  list<string> listOfFiles;
+  string directoryName;
+  directoryName = para.inPath_ + "/" + para.sourceName_;
+
+  //Make the listing of the directory
+  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
+  
+  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
+  iFileEnd = listOfFiles.end();
+  
+
+  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough...
+
+  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
+  uint_4 nSpectra=0;
+  //Loop on files
+  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
+    if (para.debuglev_>90){
+      cout << "load file <" << *iFile << ">" << endl;
+    }
+    PInPersist fin(*iFile);
+    vector<string> vec = fin.GetNameTags();
+    list<string> listOfSpectra;
+    //Keep only required PPF objects
+    std::remove_copy_if(
+			vec.begin(), vec.end(), back_inserter(listOfSpectra),
+			not1(match)
+			);
+   
+    iSpecEnd = listOfSpectra.end();
+    listOfSpectra.sort(stringCompare);
+    //Loop of spectra matrix
+    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ;  ++iSpec){
+      if (para.debuglev_>90){
+	cout << " spactra <" << *iSpec << ">" << endl;
+      }
+      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      fin.GetObject(aSpec,*iSpec);
+
+      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
+
+      nSpectra++;
+    }//eo loop on spectra in a file
+  }//eo loop on files
+  
+
+ 
+  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+  //Compute the median for each freq. and Channel
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
+    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
+      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
+      vector<r_4> tmp;
+      tmp0.FillTo(tmp);
+      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
+    }
+  }
+
+
+  //Compute the reduced version of the mean and sigma
+  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
+
+
+  sa_size_t f1320=freqTochan(1320.);
+  sa_size_t f1345=freqTochan(1345.);
+  sa_size_t f1355=freqTochan(1355.);
+  sa_size_t f1380=freqTochan(1380.);
+  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
+  if (para.debuglev_>9){
+    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
+  }
+  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    double meanMed1;
+    double sigmaMed1;
+    TVector<r_4> band1;
+    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
+    MeanSigma(band1,meanMed1,sigmaMed1);
+    double meanMed2;
+    double sigmaMed2;
+    TVector<r_4> band2;
+    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
+    MeanSigma(band2,meanMed2,sigmaMed2);
+    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
+  }  
+  meanMed.Print(cout);
+  
+  //Compute the sigma in the range 1320MHz-1380MHz
+  if (para.debuglev_>9){
+    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
+  }
+  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
+  sa_size_t redf1320=(1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
+  sa_size_t redf1380=(1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
+
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    double meanSigma;
+    double sigmaSigma;
+    TVector<r_4> band;
+    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
+    MeanSigma(band,meanSigma,sigmaSigma);
+    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
+    sigmaMed(iCh) = meanSigma;
+  }
+  sigmaMed.Print(cout);
+  
+  
+  if (para.debuglev_>9){
+    cout << "Compute medianOfSpectraNorm" << endl;
+  }
+  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
+    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
+  }
+
+  
+
+  {//Save the result
+    stringstream tmp;
+    tmp << nSpectra;
+    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
+    POutPersist fos(fileName);
+
+    string tag = "median";
+    fos << PPFNameTag(tag) << medianOfSpectra;
+
+    tag = "medianNorm";
+    fos << PPFNameTag(tag) << medianOfSpectraNorm;
+    
+
+    tag = "meanmedred";
+    fos << PPFNameTag(tag) << meanRedMtx;
+    tag = "sigmamedred";
+    fos << PPFNameTag(tag) << sigmaRedMtx;
+    tag = "cycleVsfreq";
+    
+    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
+    fos << PPFNameTag(tag) << tarr;
+  }
+}
+
+//-------------------------------------------------------
+int main(int narg, char* arg[]) {
+  
+  int rc = 0; //return code
+  string msg; //message used in Exceptions
+
+
+
+  //default value for initial parameters (see Para structure on top of the file)
+  string debuglev = "0";
+  string action = "meanDiffOnOff";
+  string inputPath = "."; 
+  string outputPath = "."; 
+  string sourceName = "Abell85";
+  string ppfFile;
+  string nSliceInFreq = "32";
+ 
+  
+  //  bool okarg=false;
+  int ka=1;
+  while (ka<narg) {
+    if (strcmp(arg[ka],"-h")==0) {
+      cout << "Usage:  -act [meanDiffOnOff]|medianDiffOnOff\n"
+	   << " -src [Abell85]\n -inPath [.]|<top_root_dir of the ppf file>\n" 
+	   << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
+	   << " -outPath [.]|<dir of the output> \n"
+	   << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
+	   << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)"
+	   << " -debug <level> "
+	   << endl;
+      return 0;
+    }
+    else if (strcmp(arg[ka],"-debug")==0) {
+      debuglev=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-act")==0) {
+      action=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-inPath")==0) {
+      inputPath=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-outPath")==0) {
+      outputPath=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-src")==0) {
+      sourceName=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-ppfFile")==0) {
+      ppfFile=arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
+      nSliceInFreq=arg[ka+1];
+      ka+=2;
+    }
+    else ka++;
+  }//eo while
+
+  para.debuglev_   = atoi(debuglev.c_str());
+  para.inPath_     = inputPath;
+  para.outPath_    = outputPath;
+  para.sourceName_ = sourceName;
+  para.ppfFile_    = ppfFile;
+  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
+
+  cout << "Dump Initial parameters ............" << endl;
+  cout << " action = " << action << "\n"
+       << " inputPath = " << para.inPath_  << "\n" 
+       << " outputPath = " <<  para.outPath_ << "\n"
+       << " sourceName = " << para.sourceName_ << "\n"
+       << " ppfFile = " <<  para.ppfFile_ << "\n"
+       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
+       << " debuglev = "  << para.debuglev_   << "\n";
+  cout << "...................................." << endl;
+
+  if ( "" == ppfFile ) {
+    cerr << "mergeRawOnOff.cc: you have forgotten ppfFile option"
+	 << endl;
+    return 999;
+  }
+
+
+  try {
+
+//     int b,e;
+//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
+//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
+
+    if ( action == "meanDiffOnOff" ) {
+      meanOnCycles();
+    } else if (action == "medianDiffOnOff") {
+      medianOnCycles();
+    } else {
+      msg = "Unknown action " + action;
+      throw msg;
+    }
+
+
+  }  catch (std::exception& sex) {
+    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
+	 << "\n msg= " << sex.what() << endl;
+    rc = 78;
+  }
+  catch ( string str ) {
+    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
+  }
+  catch (...) {
+    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
+    rc = 79; 
+  } 
+
+  return 0;
+
+}
Index: ORadio/AmasNancay/trunk/mergeRawOnOff.cc
===================================================================
--- /BAORadio/AmasNancay/trunk/mergeRawOnOff.cc	(revision 557)
+++ 	(revision )
@@ -1,546 +1,0 @@
-// Utilisation de SOPHYA pour faciliter les tests ...
-#include <regex.h>
-//#include <regexp.h>
-#include <stdio.h>
-
-#include "sopnamsp.h"
-#include "machdefs.h"
-
-#include <stdlib.h>
-#include <dirent.h>
-#include <matharr.h>
-
-// include standard c/c++
-#include <iostream>
-#include <fstream>
-#include <string>
-#include <vector>
-#include <map>
-#include <functional>
-#include <algorithm>
-#include <numeric>
-#include <list>
-#include <exception>
-
-// include sophya mesure ressource CPU/memoire ...
-#include "resusage.h"
-#include "ctimer.h"
-#include "timing.h"
-#include "timestamp.h"
-#include "strutilxx.h"
-#include "ntuple.h"
-#include "fioarr.h"
-#include "tarrinit.h"
-#include "histinit.h"
-#include "fitsioserver.h"
-#include "fiosinit.h"
-#include "ppersist.h"
-
-//-----------------------------------------------
-const sa_size_t NUMBER_OF_CHANNELS = 2;
-const sa_size_t NUMBER_OF_FREQ     = 8192;
-const sa_size_t TOTAL_NUM_CYCLES   = 500;
-const r_4    LOWER_FREQUENCY = 1250.0; //MHz
-const r_4    TOTAL_BANDWIDTH = 250.0; //MHz
-//-----------------------------------------------
-//Input parameters
-struct Param {
-  int debuglev_;      //debug
-  string inPath_;     //root directory of the input files
-  string outPath_;    //output files are located here
-  string sourceName_; //source name & subdirectory of the input files
-  string ppfFile_;    //generic name of the input files
-  int nSliceInFreq_;  //used by reduceSpectra() fnc
-} para;
-//--------------------------------------------------------------
-//Utility functions
-
-sa_size_t freqTochan(r_4 f){
-  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
-}
-//---------
-  class median_of_empty_list_exception:public std::exception{
-    virtual const char* what() const throw() {
-    return "Attempt to take the median of an empty list of numbers.  "
-      "The median of an empty list is undefined.";
-  }
-  };
-  template<class RandAccessIter>
-    double median(RandAccessIter begin, RandAccessIter end) 
-    throw(median_of_empty_list_exception){
-    if(begin == end){ throw median_of_empty_list_exception(); }
-    std::size_t size = end - begin;
-    std::size_t middleIdx = size/2;
-    RandAccessIter target = begin + middleIdx;
-    std::nth_element(begin, target, end);
-    
-    if(size % 2 != 0){ //Odd number of elements
-      return *target;
-    }else{            //Even number of elements
-      double a = *target;
-      RandAccessIter targetNeighbor= target-1;
-      std::nth_element(begin, targetNeighbor, end);
-      return (a+*targetNeighbor)/2.0;
-    }
-  }
-
-//--------------------------------------------------------------
-char *regexp (const char *string, const char *patrn, int *begin, int *end) {    
-        int i, w=0, len;                 
-        char *word = NULL;
-        regex_t rgT;
-        regmatch_t match;
-        regcomp(&rgT,patrn,REG_EXTENDED);
-        if ((regexec(&rgT,string,1,&match,0)) == 0) {
-                *begin = (int)match.rm_so;
-                *end = (int)match.rm_eo;
-                len = *end-*begin;
-                word=(char*)malloc(len+1);
-                for (i=*begin; i<*end; i++) {
-                        word[w] = string[i];
-                        w++; }
-                word[w]=0;
-        }
-        regfree(&rgT);
-        return word;
-}
-//------- 
-sa_size_t round_sa(r_4 r) {
-  return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
-}
-//-----
-string StringToLower(string strToConvert){
-  //change each element of the string to lower case
-  for(unsigned int i=0;i<strToConvert.length();i++) {
-    strToConvert[i] = tolower(strToConvert[i]);
-  }
-  return strToConvert;//return the converted string
-}
-//-----
-bool stringCompare( const string &left, const string &right ){
-   if( left.size() < right.size() )
-      return true;
-   for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
-      if( tolower( *lit ) < tolower( *rit ) )
-         return true;
-      else if( tolower( *lit ) > tolower( *rit ) )
-         return false;
-   return false;
-}
-//-----
-list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
-  list<string> theList;
-
-
-  DIR *dip;
-  struct dirent *dit;
-  string msg;  string fileName;
-  string fullFileName;
-  size_t found;
-
-  if ((dip=opendir(dir.c_str())) == NULL ) {
-    msg = "opendir failed on directory "+dir;
-    throw msg;
-  }
-  while ( (dit = readdir(dip)) != NULL ) {
-    fileName = dit->d_name;
-    found=fileName.find(filePettern);
-    if (found != string::npos) {
-      fullFileName = dir + "/";
-      fullFileName += fileName;
-      theList.push_back(fullFileName);
-    }
-  }//eo while
-  if (closedir(dip) == -1) {
-    msg = "closedir failed on directory "+dir;
-    throw msg;
-  }
-  
-  theList.sort(stringCompare);
-
-  return theList;
-
-}
-//
-class  StringMatch : public unary_function<string,bool> {
-public:
-  StringMatch(const string& pattern): pattern_(pattern){}
-  bool operator()(const string& aStr) const {
-
-
-    int b,e;
-    regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
-
-//     cout << "investigate " << aStr << " to find " << pattern_ 
-// 	 << "[" <<b<<","<<e<<"]" 
-// 	 << endl;
-
-    
-    if (b != 0) return false;
-    if (e != (int)aStr.size()) return false;
-    return true;
-
-  }
-private:
-  string pattern_;
-};
-//-------------------------------------------------------
-//Rebin in frequence + compute mean and sigma
-void reduceSpectra(const TMatrix<r_4>& specMtxInPut, 
-		    TMatrix<r_4>& meanMtx, 
-		    TMatrix<r_4>& sigmaMtx) {
-  sa_size_t nSliceFreq = para.nSliceInFreq_;
-  sa_size_t deltaFreq =  NUMBER_OF_FREQ/nSliceFreq;
-  for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
-    sa_size_t freqLow= iSlice*deltaFreq;
-    sa_size_t freqHigh= freqLow + deltaFreq -1;
-    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
-      TVector<r_4> reducedRow;
-      reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
-      double mean; 
-      double sigma;
-      MeanSigma(reducedRow,mean,sigma);
-      meanMtx(iCh,iSlice) = mean;
-      sigmaMtx(iCh,iSlice) = sigma;
-    }//eo loop on channels
-  }//eo loop on slices
-}
-//-------------------------------------------------------
-//Compute the mean of Raw spectra and also the mean/sigma of rebinned spectra 
-//Used like:
-//
-void meanOnCycles() throw(string) {
-  list<string> listOfFiles;
-  string directoryName;
-  directoryName = para.inPath_ + "/" + para.sourceName_;
-
-  //Make the listing of the directory
-  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
-  
-  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
-  iFileEnd = listOfFiles.end();
-  
-  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
-  TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-  uint_4 nSpectra=0;
-  //Loop on files
-  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
-    if (para.debuglev_>90){
-      cout << "load file <" << *iFile << ">" << endl;
-    }
-    PInPersist fin(*iFile);
-    vector<string> vec = fin.GetNameTags();
-    list<string> listOfSpectra;
-    //Keep only required PPF objects
-    std::remove_copy_if(
-			vec.begin(), vec.end(), back_inserter(listOfSpectra),
-			not1(match)
-			);
-   
-    iSpecEnd = listOfSpectra.end();
-    listOfSpectra.sort(stringCompare);
-    //Loop of spectra matrix
-    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
-      if (para.debuglev_>90){
-	cout << " spactra <" << *iSpec << ">" << endl;
-      }
-      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-      fin.GetObject(aSpec,*iSpec);
-      //How to see if the GetObject is ok?? Ask Reza
-      nSpectra++;
-      meanOfSpectra+=aSpec;
-    }//eo loop on spectra in a file
-  }//eo loop on files
-  
-  //Normalisation
-  if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
-
-  //Compute the reduced version of the mean and sigma
-  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
-  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
-  reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
-
-  {//Save the result
-    stringstream tmp;
-    tmp << nSpectra;
-    string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
-    cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
-    POutPersist fos(fileName);
-
-    string tag = "mean";
-    fos << PPFNameTag(tag) << meanOfSpectra;
-    tag = "meanred";
-    fos << PPFNameTag(tag) << meanRedMtx;
-    tag = "sigmared";
-    fos << PPFNameTag(tag) << sigmaRedMtx;
-  }
-}
-//-------------------------------------------------------
-//Compute the median of Raw spectra and also the mean/sigma of rebinned spectra 
-//Used like:
-//
-void medianOnCycles() throw(string) {
-  list<string> listOfFiles;
-  string directoryName;
-  directoryName = para.inPath_ + "/" + para.sourceName_;
-
-  //Make the listing of the directory
-  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
-  
-  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
-  iFileEnd = listOfFiles.end();
-  
-
-  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough...
-
-  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
-  uint_4 nSpectra=0;
-  //Loop on files
-  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
-    if (para.debuglev_>90){
-      cout << "load file <" << *iFile << ">" << endl;
-    }
-    PInPersist fin(*iFile);
-    vector<string> vec = fin.GetNameTags();
-    list<string> listOfSpectra;
-    //Keep only required PPF objects
-    std::remove_copy_if(
-			vec.begin(), vec.end(), back_inserter(listOfSpectra),
-			not1(match)
-			);
-   
-    iSpecEnd = listOfSpectra.end();
-    listOfSpectra.sort(stringCompare);
-    //Loop of spectra matrix
-    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ;  ++iSpec){
-      if (para.debuglev_>90){
-	cout << " spactra <" << *iSpec << ">" << endl;
-      }
-      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-      fin.GetObject(aSpec,*iSpec);
-
-      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
-
-      nSpectra++;
-    }//eo loop on spectra in a file
-  }//eo loop on files
-  
-
- 
-  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-  //Compute the median for each freq. and Channel
-  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
-    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
-      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
-      vector<r_4> tmp;
-      tmp0.FillTo(tmp);
-      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
-    }
-  }
-
-
-  //Compute the reduced version of the mean and sigma
-  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
-  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
-  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
-
-
-  sa_size_t f1320=freqTochan(1320.);
-  sa_size_t f1345=freqTochan(1345.);
-  sa_size_t f1355=freqTochan(1355.);
-  sa_size_t f1380=freqTochan(1380.);
-  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
-  if (para.debuglev_>9){
-    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
-  }
-  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
-  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
-    double meanMed1;
-    double sigmaMed1;
-    TVector<r_4> band1;
-    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
-    MeanSigma(band1,meanMed1,sigmaMed1);
-    double meanMed2;
-    double sigmaMed2;
-    TVector<r_4> band2;
-    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
-    MeanSigma(band2,meanMed2,sigmaMed2);
-    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
-  }  
-  meanMed.Print(cout);
-  
-  //Compute the sigma in the range 1320MHz-1380MHz
-  if (para.debuglev_>9){
-    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
-  }
-  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
-  sa_size_t redf1320=(1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
-  sa_size_t redf1380=(1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
-
-  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
-    double meanSigma;
-    double sigmaSigma;
-    TVector<r_4> band;
-    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
-    MeanSigma(band,meanSigma,sigmaSigma);
-    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
-    sigmaMed(iCh) = meanSigma;
-  }
-  sigmaMed.Print(cout);
-  
-  
-  if (para.debuglev_>9){
-    cout << "Compute medianOfSpectraNorm" << endl;
-  }
-  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
-  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
-    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
-    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
-  }
-
-  
-
-  {//Save the result
-    stringstream tmp;
-    tmp << nSpectra;
-    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
-    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
-    POutPersist fos(fileName);
-
-    string tag = "median";
-    fos << PPFNameTag(tag) << medianOfSpectra;
-
-    tag = "medianNorm";
-    fos << PPFNameTag(tag) << medianOfSpectraNorm;
-    
-
-    tag = "meanmedred";
-    fos << PPFNameTag(tag) << meanRedMtx;
-    tag = "sigmamedred";
-    fos << PPFNameTag(tag) << sigmaRedMtx;
-    tag = "cycleVsfreq";
-    
-    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
-    fos << PPFNameTag(tag) << tarr;
-  }
-}
-
-//-------------------------------------------------------
-int main(int narg, char* arg[]) {
-  
-  int rc = 0; //return code
-  string msg; //message used in Exceptions
-
-
-
-  //default value for initial parameters (see Para structure on top of the file)
-  string debuglev = "0";
-  string action = "meanDiffOnOff";
-  string inputPath = "."; 
-  string outputPath = "."; 
-  string sourceName = "Abell85";
-  string ppfFile;
-  string nSliceInFreq = "32";
- 
-  
-  //  bool okarg=false;
-  int ka=1;
-  while (ka<narg) {
-    if (strcmp(arg[ka],"-h")==0) {
-      cout << "Usage:  -act [meanDiffOnOff]|medianDiffOnOff\n"
-	   << " -src [Abell85]\n -inPath [.]|<top_root_dir of the ppf file>\n" 
-	   << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
-	   << " -outPath [.]|<dir of the output> \n"
-	   << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
-	   << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)"
-	   << " -debug <level> "
-	   << endl;
-      return 0;
-    }
-    else if (strcmp(arg[ka],"-debug")==0) {
-      debuglev=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-act")==0) {
-      action=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-inPath")==0) {
-      inputPath=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-outPath")==0) {
-      outputPath=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-src")==0) {
-      sourceName=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-ppfFile")==0) {
-      ppfFile=arg[ka+1];
-      ka+=2;
-    }
-    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
-      nSliceInFreq=arg[ka+1];
-      ka+=2;
-    }
-    else ka++;
-  }//eo while
-
-  para.debuglev_   = atoi(debuglev.c_str());
-  para.inPath_     = inputPath;
-  para.outPath_    = outputPath;
-  para.sourceName_ = sourceName;
-  para.ppfFile_    = ppfFile;
-  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
-
-  cout << "Dump Initial parameters ............" << endl;
-  cout << " action = " << action << "\n"
-       << " inputPath = " << para.inPath_  << "\n" 
-       << " outputPath = " <<  para.outPath_ << "\n"
-       << " sourceName = " << para.sourceName_ << "\n"
-       << " ppfFile = " <<  para.ppfFile_ << "\n"
-       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
-       << " debuglev = "  << para.debuglev_   << "\n";
-  cout << "...................................." << endl;
-
-  if ( "" == ppfFile ) {
-    cerr << "mergeRawOnOff.cc: you have forgotten ppfFile option"
-	 << endl;
-    return 999;
-  }
-
-
-  try {
-
-//     int b,e;
-//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
-//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
-
-    if ( action == "meanDiffOnOff" ) {
-      meanOnCycles();
-    } else if (action == "medianDiffOnOff") {
-      medianOnCycles();
-    } else {
-      msg = "Unknown action " + action;
-      throw msg;
-    }
-
-
-  }  catch (std::exception& sex) {
-    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
-	 << "\n msg= " << sex.what() << endl;
-    rc = 78;
-  }
-  catch ( string str ) {
-    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
-  }
-  catch (...) {
-    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
-    rc = 79; 
-  } 
-
-  return 0;
-
-}
