source: BAORadio/AmasNancay/trunk/mergeAnaFiles.cc @ 606

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"
#include <dirent.h>
#include <matharr.h>

// include standard c/c++
#include <regex.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits>
#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"

//-----------------------------------------------
//Usage
//
//./Objs/mergeAnaFiles -act meanRawDiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src NGC4383 -ppfFile dataRaw -debug 1000 -mxcycles 10
//./Objs/mergeAnaFiles -src Abell85 -act meanCalibBAODiffOnOff -mxcycles 500 -calibfreq 1410 -inPath /sps/baoradio/AmasNancay/JEC/ -ppfFile dataRaw -debug 1 >& mergeAna-500.log
//
//
//-----------------------------------------------
const sa_size_t NUMBER_OF_CHANNELS = 2;
const sa_size_t NUMBER_OF_FREQ     = 8192;
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
  string calibFreq_;  //freq. value used for calibration
  r_4 rcalibFreq_;    //float version
  string calibBandFreq_;  //band of freq. used for calibration
  r_4 rcalibBandFreq_;    //float version
  int maxNumberCycles_;//maximum number of cycles to be treated
} para;
//--------------------------------------------------------------
//Utility functions

sa_size_t freqToChan(r_4 f){
  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
}
//--------
//COMPUTE the mean value on a freq. range for all channels
//--------
void meanInRange(const TMatrix<r_4> mtx, 
                 sa_size_t chLow,
                 sa_size_t chHigh, 
                 TVector<r_4>& vec){
  
  sa_size_t nr = mtx.NRows();
  
  for (sa_size_t ir=0; ir<nr; ir++){
    TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
    double mean, sigma;
    MeanSigma(tmp,mean,sigma);
    vec(ir) = mean;
  }
}
//---------
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;
  }
}

//-------------
//JEC 25/10/11 Perform a median filtering with a sliding window of "halfwidth" half width
//             It takes care of the edges and is based on the median function (above) 
void medianFiltering(const TMatrix<r_4> mtx, 
                     sa_size_t halfwidth,
                     TMatrix<r_4>& vec) {
  
  sa_size_t nr = mtx.NRows();
  sa_size_t nc = mtx.NCols();
  sa_size_t chMin = 0;
  sa_size_t chMax = nc-1;
  
  for (sa_size_t ir=0; ir<nr; ir++){
    for (sa_size_t ic=0; ic<nc; ic++) {
      sa_size_t chLow = ic-halfwidth;
      chLow = (chLow >= chMin) ? chLow : chMin;
      sa_size_t chHigh = ic+halfwidth;
      chHigh = ( chHigh <= chMax ) ? chHigh : chMax;
      TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
      vector<r_4> val;
      tmp.FillTo(val);
      vec(ir,ic) = median(val.begin(),val.end());
    }
  }
}
//-------------
void split(const string& str, const string& delimiters , vector<string>& tokens) {
    // Skip delimiters at beginning.
    string::size_type lastPos = str.find_first_not_of(delimiters, 0);
    // Find first "non-delimiter".
    string::size_type pos     = str.find_first_of(delimiters, lastPos);

    while (string::npos != pos || string::npos != lastPos)
    {
        // Found a token, add it to the vector.
        tokens.push_back(str.substr(lastPos, pos - lastPos));
        // Skip delimiters.  Note the "not_of"
        lastPos = str.find_first_not_of(delimiters, pos);
        // Find next "non-delimiter"
        pos = str.find_first_of(delimiters, lastPos);
    }
}
//--------------------------------------------------------------
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 Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra 
//
void meanCalibBAODiffOnOffCycles() 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();

  //mean ON-OFF over the list of cycles
  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0 
  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0 
  TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
  TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
  static const int NINFO=25;
  char* onffTupleName[NINFO]={"cycle"
                              ,"onoffRaw0","onoffRaw1"
                              ,"onoffRun0","onoffRun1"
                              ,"onoffCycle0","onoffCycle1"
                              ,"onoffRaw01420","onoffRaw11420"
                              ,"onoffRun01420","onoffRun11420"
                              ,"onoffCycle01420","onoffCycle11420"
                              ,"onoffRaw01420side","onoffRaw11420side"
                              ,"onoffRun01420side","onoffRun11420side"
                              ,"onoffCycle01420side","onoffCycle11420side"
                              ,"onoffRaw0f14001420","onoffRaw1f14001420"
                              ,"offRaw0f14001420","offRaw1f14001420"
                              ,"onRaw0f14001420","onRaw1f14001420"
  };
  NTuple onoffevolution(NINFO,onffTupleName);
  r_4 xnt[NINFO];

  //Lower and Higher freq. arround the Calibration Freq. bin to perform mean follow up
  sa_size_t chCalibLow  = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
  sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
  sa_size_t ch1420High = freqToChan(1420.4+0.2);

  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
  sa_size_t ch1420aLow  = freqToChan(1418);
  sa_size_t ch1420aHigh = freqToChan(1419);
  sa_size_t ch1420bLow  = freqToChan(1422);
  sa_size_t ch1420bHigh = freqToChan(1423);
  
  //1400-1420Mhz
  sa_size_t ch1400 = freqToChan(1400);
  sa_size_t ch1405 = freqToChan(1400);
  sa_size_t ch1415 = freqToChan(1415);
  sa_size_t ch1420 = freqToChan(1420);

  if (para.debuglev_>0){
    cout << "freq. band for follow up [" <<  chCalibLow << ", "<< chCalibHigh << "]" << endl;
    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
  }
  
  //Loop on files/run

  int totalNumberCycles=0; //total number of cycles for normalisation
  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
    if (para.debuglev_>90){
      cout << "load file <" << *iFile << ">" << endl;
    }

    vector<string> tokens;
    split(*iFile,"_",tokens);
    string dateOfRun = tokens[1];
    if (para.debuglev_>90){
      cout << "date <" << dateOfRun << ">" << endl;
    }
    vector<string> tokens2;
    split(tokens[2],".",tokens2);
    string srcLower = tokens2[0];


    
    PInPersist fin(*iFile);
    vector<string> vec = fin.GetNameTags();

    vector<string> modeList;
    modeList.push_back("On");
    modeList.push_back("Off");
    vector<string>::const_iterator iMode;
    
    map<string, list<int> > cycleModeCollect;
    
    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
      list<string> listOfSpectra;
      //Keep only required PPF objects
      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
      std::remove_copy_if(
                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
                          not1(StringMatch(matchstr))
                          );
      
      listOfSpectra.sort(stringCompare);
      iSpecEnd = listOfSpectra.end();
      
      matchstr = "[0-9]+";
      //Loop of spectra matrix
      list<int> listOfCycles;
      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
        int b,e;
        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
        if (para.debuglev_>90){
          cout << " spactra <" << *iSpec << ">" << endl;
          cout << " cycle " << iSpec->substr(b) << endl;
        }
        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
      }//end loop spectra
      cycleModeCollect[*iMode] = listOfCycles;
    }//end of mode    

    //Take the Intersection of the list Of cycles in mode Off and On    
    list<int> commonCycles;
    set_intersection(cycleModeCollect["On"].begin(),
                     cycleModeCollect["On"].end(),
                     cycleModeCollect["Off"].begin(),
                     cycleModeCollect["Off"].end(),
                     back_inserter(commonCycles)
                     );
    
    if (para.debuglev_>90){
      cout << "Liste of cycles common to On & Off: <";
      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
        cout << *i << " ";
      }
      cout << ">" << endl;
    }
    
    //
    //Load BAO Calibration factors "per Cycle and Channels"
    //Compute the mean per Cycle to 
    //  fill factors "per Run and Channels" with the same cycle list
    //
    //
    //TODO improve the code....

    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
    
    string calibFileName;
    ifstream ifs;
    sa_size_t nr,nc; //values read 

    //OFF Cycle per Channel
    calibFileName = directoryName + "/" 
      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
    ifs.open(calibFileName.c_str());
    if ( ! ifs.is_open() ) {

      throw calibFileName + " cannot be opened...";
    }   
    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
    if(para.debuglev_>9){
      cout << "(nr,nc): "<< nr << "," << nc << endl;
      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
      cout << endl;
    }

    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
    {//Compute the mean
      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
      double mean,sigma;
      MeanSigma(coef,mean,sigma);
      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
    }
    if(para.debuglev_>9){
      cout << "Fill calib. with mean value " << endl; 
      calibBAOfactors_Off_Run_Ch0.Print(cout);
      cout << endl;
    }
    ifs.close();

    //
    calibFileName = directoryName + "/" 
      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
    ifs.open(calibFileName.c_str());
    if ( ! ifs.is_open() ) {

      throw calibFileName + " cannot be opened...";
    }   
    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
    if(para.debuglev_>9){
      cout << "(nr,nc): "<< nr << "," << nc << endl;
      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
      cout << endl;
    }
    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
    {//Compute the mean
      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
      double mean,sigma;
      MeanSigma(coef,mean,sigma);
      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
    }
    if(para.debuglev_>9){
      cout << "Fill calib. with mean value " << endl; 
      calibBAOfactors_Off_Run_Ch1.Print(cout);
      cout << endl;
    }
    ifs.close();

    //ON Cycle per Channel
    calibFileName = directoryName + "/" 
      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
    ifs.open(calibFileName.c_str());
    if ( ! ifs.is_open() ) {

      throw calibFileName + " cannot be opened...";
    }   
    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
    if(para.debuglev_>9){
      cout << "(nr,nc): "<< nr << "," << nc << endl;
      calibBAOfactors_On_Cycle_Ch0.Print(cout);
      cout << endl;      
    }

    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
    {//Compute the mean
      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
      double mean,sigma;
      MeanSigma(coef,mean,sigma);
      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
    }
    if(para.debuglev_>9){
      cout << "Fill calib. with mean value " << endl; 
      calibBAOfactors_On_Run_Ch0.Print(cout);
      cout << endl;
    }
    ifs.close();

    
    calibFileName = directoryName + "/" 
      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
    ifs.open(calibFileName.c_str());
    if ( ! ifs.is_open() ) {
      throw calibFileName + " cannot be opened...";
    }   
    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
    if(para.debuglev_>9){
      cout << "(nr,nc): "<< nr << "," << nc << endl;
      calibBAOfactors_On_Cycle_Ch1.Print(cout);
      cout << endl;
    }
    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
    {//Compute the mean
      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
      double mean,sigma;
      MeanSigma(coef,mean,sigma);
      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
    }
    if(para.debuglev_>9){
      cout << "Fill calib. with mean value " << endl; 
      calibBAOfactors_On_Run_Ch1.Print(cout);
      cout << endl;
    }

    ifs.close();
    
    //link <cycle> - <calibration coefficient>
    //We cannot rely on identical cycle list of the OFF and ON calibration
    map<int,r_4> calibBAO_Off_Run_Ch0;
    map<int,r_4> calibBAO_Off_Run_Ch1;
    map<int,r_4> calibBAO_On_Run_Ch0;
    map<int,r_4> calibBAO_On_Run_Ch1;

    map<int,r_4> calibBAO_Off_Cycle_Ch0;
    map<int,r_4> calibBAO_Off_Cycle_Ch1;
    map<int,r_4> calibBAO_On_Cycle_Ch0;
    map<int,r_4> calibBAO_On_Cycle_Ch1;

    //per Run based BAO coefficients
    nr = calibBAOfactors_Off_Run_Ch0.NRows();
    for (sa_size_t ir=0; ir<nr; ++ir){
//       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
//         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;

      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
        = calibBAOfactors_Off_Run_Ch0(ir,1);
      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
        = calibBAOfactors_Off_Run_Ch1(ir,1);
      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
    }//eo loop on coef Off
    
    nr = calibBAOfactors_On_Run_Ch0.NRows();
    for (sa_size_t ir=0; ir<nr; ++ir){
      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
        = calibBAOfactors_On_Run_Ch0(ir,1);
      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
        = calibBAOfactors_On_Cycle_Ch0(ir,1);
      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
        = calibBAOfactors_On_Run_Ch1(ir,1);
      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
        = calibBAOfactors_On_Cycle_Ch1(ir,1);
    }//eo loop on coef On
      
    //Loop on cyles
    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){

      //Look if the cycle has been calibrated...
      bool isCycleCalibrated = 
        ( calibBAO_On_Run_Ch0.count(*ic)    *
          calibBAO_On_Run_Ch1.count(*ic)    *
          calibBAO_Off_Run_Ch0.count(*ic)   *
          calibBAO_Off_Run_Ch1.count(*ic)   *
          calibBAO_On_Cycle_Ch0.count(*ic)  *
          calibBAO_On_Cycle_Ch1.count(*ic)  *
          calibBAO_Off_Cycle_Ch0.count(*ic) *
          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;

      if(para.debuglev_>9) {
        cout << "Calibration coefficients for cycle "<<*ic << endl; 
        if (isCycleCalibrated) {
          cout << "Cycle calibrated " << endl;
          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
        } else {
          cout << "Cycle NOT calibrated " << endl;
        }
      }//debug


      if ( ! isCycleCalibrated ) continue;
      
      string ppftag;
      //load ON phase
      stringstream cycle;
      cycle << *ic;
      
      ppftag = "specRawOn"+cycle.str();
      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      fin.GetObject(aSpecOn,ppftag);

      TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
      aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
      aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];

      TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
      aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
      aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
      
      //Load OFF phase
      ppftag = "specRawOff"+cycle.str();
      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      fin.GetObject(aSpecOff,ppftag);

      TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
      aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
      aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];

      TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
      aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
      aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];


      //Perform the difference ON-OFF with the different calibration options
      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
      meanDiffONOFF_noCalib += diffOnOff_noCalib;
      
      //JEC 29/10/11 add ON-OFF/OFF
      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
      
      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;

      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;

      TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;

      //
      totalNumberCycles++;
      //Fill NTuple
      xnt[0] = totalNumberCycles;
 
      //Follow up arround the Calibration Freq.
      TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
      xnt[1] = meanInRange_CalibFreq_noCalib(0);
      xnt[2] = meanInRange_CalibFreq_noCalib(1);

      TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
      xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
      xnt[4] = meanInRange_CalibFreq_perRunCalib(1);

      TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
      xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
      xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);


      //Follow up arround the 1420.4MHz Freq.
      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
      xnt[7] = meanInRange_1420Freq_noCalib(0);
      xnt[8] = meanInRange_1420Freq_noCalib(1);

      TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
      xnt[9] = meanInRange_1420Freq_perRunCalib(0);
      xnt[10] = meanInRange_1420Freq_perRunCalib(1);

      TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
      xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
      xnt[12] = meanInRange_1420Freq_perCycleCalib(1);


      //Follow up below the 1420.4MHz Freq.
      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);

      xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
      xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;


      TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
      TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);

      xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
      xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;


      TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
      TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);

      xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
      xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;


      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);

      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
      aSpecOffovOff.Div(aSpecOffFitltered);

      TVector<r_4> meanInRange_1400a1420Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(aSpecOffovOff,ch1400,ch1420,meanInRange_1400a1420Freq_OFF_noCalib);

      xnt[21] = meanInRange_1400a1420Freq_OFF_noCalib(0);
      xnt[22] = meanInRange_1400a1420Freq_OFF_noCalib(1);

      //      TMatrix<r_4> aSpecOnFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      //      medianFiltering(aSpecOn,halfWidth,aSpecOnFitltered);
      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
      aSpecOnovOff.Div(aSpecOffFitltered);

      TVector<r_4> meanInRange_1400a1420Freq_ON_noCalib(NUMBER_OF_CHANNELS);
      meanInRange(aSpecOnovOff,ch1400,ch1420,meanInRange_1400a1420Freq_ON_noCalib);

      xnt[23] = meanInRange_1400a1420Freq_ON_noCalib(0);
      xnt[24] = meanInRange_1400a1420Freq_ON_noCalib(1);
      
      //store infos to Ntuple
      onoffevolution.Fill(xnt);

      //Quit if enough
      if (totalNumberCycles >= para.maxNumberCycles_) break;    

    }//eo loop on cycles
    if (totalNumberCycles >= para.maxNumberCycles_) break;          
  
  }//eo loop on spectra in a file
  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
  //Normalisation
  if(totalNumberCycles > 0){
    //JEC 29/10 add ON-OFF/OFF
    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
    meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
  }  
  
  //Compute the reduced version of the mean and sigma
  TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);

  TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);

  TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
  reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);

  {//save the results
    stringstream tmp;
    tmp << totalNumberCycles;
    string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";

    POutPersist fos(fileName);
    cout << "Save output in " << fileName << endl;

    string tag = "meanNoCalib";
    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
    
    //JEC 29/10/11
    tag = "meanOvOffNoCalib";
    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;

    tag = "meanPerRunCalib";
    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
    tag = "meanPerCycleCalib";
    fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;

    tag = "redmeanNoCalib";
    fos << PPFNameTag(tag) << meanRedMtx_noCalib;
    tag = "redsigmaNoCalib";
    fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;

    tag = "redmeanPerRunCalib";
    fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
    tag = "redsigmaPerRunCalib";
    fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;

    tag = "redmeanPerCycleCalib";
    fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
    tag = "redsigmaPerCycleCalib";
    fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
    
    tag = "onoffevol";
    fos << PPFNameTag(tag) << onoffevolution;    
  }//end of save
}
//JEC 14/11/11 New meanRawDiffOnOffCycles START
//-------------------------------------------------------
//Compute the mean of Diff ON-OFF/OFF from Raw spectra
//Used like:
//
void meanRawDiffOnOffCycles() 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();

  //mean ON-OFF over the list of cycles
  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF

  int totalNumberCycles=0; //total number of cycles

  //Loop on files
  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
    if (para.debuglev_>90){
    }
    PInPersist fin(*iFile);
    vector<string> vec = fin.GetNameTags();

    vector<string> modeList;
    modeList.push_back("On");
    modeList.push_back("Off");
    vector<string>::const_iterator iMode;
    
    map<string, list<int> > cycleModeCollect;
    
    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
      list<string> listOfSpectra;
      //Keep only required PPF objects
      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
      std::remove_copy_if(
                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
                          not1(StringMatch(matchstr))
                          );
      
      listOfSpectra.sort(stringCompare);
      iSpecEnd = listOfSpectra.end();
      
      matchstr = "[0-9]+";
      //Loop of spectra matrix
      list<int> listOfCycles;
      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
        int b,e;
        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
        if (para.debuglev_>90){
          cout << " spactra <" << *iSpec << ">" << endl;
          cout << " cycle " << iSpec->substr(b) << endl;
        }
        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
      }//end loop spectra
      cycleModeCollect[*iMode] = listOfCycles;
    }//end of mode    

    //Take the Intersection of the list Of cycles in mode Off and On    
    list<int> commonCycles;
    set_intersection(cycleModeCollect["On"].begin(),
                     cycleModeCollect["On"].end(),
                     cycleModeCollect["Off"].begin(),
                     cycleModeCollect["Off"].end(),
                     back_inserter(commonCycles)
                     );
    
    if (para.debuglev_>90){
      cout << "Liste of cycles common to On & Off: <";
      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
        cout << *i << " ";
      }
      cout << ">" << endl;
    }

    //Loop on cyles
    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
      
      string ppftag;
      //load ON phase
      stringstream cycle;
      cycle << *ic;
      ppftag = "specRawOn"+cycle.str();
      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      fin.GetObject(aSpecOn,ppftag);

      //Load OFF phase
      ppftag = "specRawOff"+cycle.str();
      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      fin.GetObject(aSpecOff,ppftag);
      
      //Perform the difference ON-OFF
      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
      meanDiffONOFF_noCalib += diffOnOff_noCalib;
      
      //JEC 29/10/11 add ON-OFF/OFF
      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
      
      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;


      //JEC 15/11/11 add ON/OFF and OFF/OFF
      TMatrix<r_4> onOvOff(aSpecOn,false);
      onOvOff.Div(aSpecOffFitltered);
      meanONovOFF_noCalib += onOvOff;
      
      TMatrix<r_4> offOvOff(aSpecOff,false);
      offOvOff.Div(aSpecOffFitltered);
      meanOFFovOFF_noCalib += offOvOff;

      totalNumberCycles++;

      //Quit if enough
      if (totalNumberCycles >= para.maxNumberCycles_) break;    
    }//end of cycles

    if (totalNumberCycles >= para.maxNumberCycles_) break;          

  }//end files
  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
  //Normalisation
  if(totalNumberCycles > 0){
    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
    meanONovOFF_noCalib         /= (r_4)totalNumberCycles;
    meanOFFovOFF_noCalib        /= (r_4)totalNumberCycles;
  }  

  {//save results
    stringstream tmp;
    tmp << totalNumberCycles;
    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";

    POutPersist fos(fileName);
    cout << "Save output in " << fileName << endl;

    string tag = "meanNoCalib";
    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
    
    tag = "meanOvOffNoCalib";
    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
    
    tag = "meanOnovOffNoCalib";
    fos << PPFNameTag(tag) << meanONovOFF_noCalib;

    tag = "meanOffovOffNoCalib";
    fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;

  }//end save
}
//JEC 14/11/11 New meanRawDiffOnOffCycles END
//-------------------------------------------------------
//JEC 14/11/11 Obsolete START
//-------------------------------------------------------
//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra 
//Used like:
//
// void meanRawDiffOnOffCycles() 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)
//                      );
   
//     listOfSpectra.sort(stringCompare);
//     iSpecEnd = listOfSpectra.end();
//     //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;
//   }
// }
//JEC 14/11/11 Obsolete END
//-------------------------------------------------------
//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra 
//Used like:
//
void medianRawDiffOnOffCycles() 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,para.maxNumberCycles_); //para.maxNumberCycles_ 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)
                        );
   
    listOfSpectra.sort(stringCompare);
    iSpecEnd = listOfSpectra.end();
    //Loop of spectra matrix
    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ;  ++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);
  cout << endl;

  
  //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=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
  sa_size_t redf1380=(sa_size_t)((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);
  cout << endl;

  
  
  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";
  string typeOfCalib="perRun";
  string calibFreq = "1346";
  string calibBandFreq="6.25";
  string mxcycles;

  //  bool okarg=false;
  int ka=1;
  while (ka<narg) {
    if (strcmp(arg[ka],"-h")==0) {
      cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
           << " -mxcycles <number> (max. number of cycles to be treated)\n"
           << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
           << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\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)\n"
           << " -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],"-calibfreq")==0) {
      calibFreq=arg[ka+1];
      ka+=2;
    }    
    else if (strcmp(arg[ka],"-calibbandfreq")==0) {
      calibBandFreq=arg[ka+1];
      ka+=2;
    }    
    else if (strcmp(arg[ka],"-mxcycles")==0) {
      mxcycles=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());
  para.calibFreq_   = calibFreq;
  para.calibBandFreq_ = calibBandFreq;
  para.rcalibFreq_   = atof(calibFreq.c_str());
  para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
  if (mxcycles != "") {
    para.maxNumberCycles_ = atoi(mxcycles.c_str());
  } else {
    para.maxNumberCycles_ = std::numeric_limits<int>::max();
  }

  cout << "Dump Initial parameters ............" << endl;
  cout << " action = " << action << "\n"
       << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
       << " inputPath = " << para.inPath_  << "\n" 
       << " outputPath = " <<  para.outPath_ << "\n"
       << " sourceName = " << para.sourceName_ << "\n"
       << " ppfFile = " <<  para.ppfFile_ << "\n"
       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
       << " calibFreq = " <<  para.calibFreq_ << "\n"
       << " calibBandFreq = " <<  para.calibBandFreq_ << "\n"
       << " debuglev = "  << para.debuglev_   << "\n";
  cout << "...................................." << endl;

  if ( "" == ppfFile ) {
    cerr << "mergeAnaFiles.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 == "meanRawDiffOnOff" ) {
      meanRawDiffOnOffCycles();
    } else if (action == "medianRawDiffOnOff") {
      medianRawDiffOnOffCycles();
    } else if (action == "meanCalibBAODiffOnOff") {
      meanCalibBAODiffOnOffCycles();
    } 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;

}