source: BAORadio/AmasNancay/analyse.cc@ 523

//
// Analyse of Amas@Nancay runs by J.E Campagne (LAL)
// Version 0: 1/6/2011
//-----------------------------

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

#include <stdlib.h>
#include <dirent.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 r_4    LOWER_FREQUENCY = 1250.0; //MHz
const r_4    TOTAL_BANDWIDTH = 250.0; //MHz

//decalration of non class members functions
extern "C" {
  int Usage(bool);
}


//----------------------------------------------------------
//Utility fonctions
// Function for deleting pointers in map.
template<class A, class B>
struct DeleteMapFntor
{
    // Overloaded () operator.
    // This will be called by for_each() function.
    bool operator()(pair<A,B> x) const
    {
        // Assuming the second item of map is to be
        // deleted. Change as you wish.
        delete x.second;
        return true;
    }
};
//-----
bool compare(const pair<string,r_4>& i, const pair<string,r_4>& j) {
  return i.second < j.second;
}
//-----
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
}
//-----
//JEC 22/9/11 comparison, not case sensitive to sort File liste START
bool stringCompare( const string &left, const string &right ){
   if( left.size() < right.size() )
      return true;
   return false;
   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;
}//JEC 22/9/11 comparison, not case sensitive to sort File liste END
//-----
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;
  }
  
  //JEC 22/9/11 START
  theList.sort(stringCompare);
  //JEC 22/9/11 END

  return theList;
}


//Declaration of local classes
//----------------------------------------------
//Process Interface
class IProcess {
public:
  IProcess() {}
  virtual ~IProcess() {}
  virtual int processCmd() throw(string) =0;
};
//------------
//Common Process
class ProcessBase : public IProcess {
public:
  ProcessBase();
  virtual ~ProcessBase();
  void SetInputPath(const string& inputPath) {inputPath_ = inputPath;}
  void SetOutputPath(const string& outputPath) {outputPath_ = outputPath;}
  void SetSourceName(const string& sourceName) {sourceName_ = sourceName;}
  void SetDateOfRun(const string& dateOfRun) {dateOfRun_ = dateOfRun;}
  void SetSpectraDirectory(const string& spectraDirectory) {spectraDirectory_ = spectraDirectory;}
  void SetTypeOfFile(const string& typeOfFile) { typeOfFile_ = typeOfFile; }
  void SetNumCycle(const string& numcycle) {numcycle_ = numcycle; }
  void SetScaFileName(const string& scaFileName) { scaFileName_ =scaFileName; }

  void SetDebugLevel(const string& debuglev) {
    debuglev_ = atoi(debuglev.c_str());
  }

  virtual int processCmd() throw(string);
  
protected:
  string inputPath_;
  string outputPath_;
  string sourceName_;
  string dateOfRun_;
  string spectraDirectory_;
  string typeOfFile_;

  string numcycle_; //cycle numbers format="first,last"
  sa_size_t ifirstCycle_;
  sa_size_t ilastCycle_;


  uint_4 debuglev_; 
  string scaFileName_;
  NTuple* scaTuple_;
  map<sa_size_t,sa_size_t> idCycleInTuple_;
};
ProcessBase::ProcessBase() {
  scaTuple_ = 0;
}
ProcessBase::~ProcessBase() {
  if (scaTuple_) delete scaTuple_;
  scaTuple_ = 0;
}
//------------
//Process ON/OFF data
//------------
class ProcessONOFFData : public ProcessBase {
protected:
  string  freqBAOCalibration_;//string MHz
public:
  ProcessONOFFData(){}
  virtual ~ProcessONOFFData(){}

  void SetFreqBAOCalibration(const string& freqBAOCalibration) { 
    freqBAOCalibration_ = freqBAOCalibration; 
  }
  
  virtual int processCmd() throw(string);
};
//------------
//Process Gain
//------------
class ProcessGain : public ProcessBase {
protected:
  string mode_; //mode of data taken for gain computation On || Off
public:
  ProcessGain(){}
  virtual ~ProcessGain(){}

  void SetMode(const string& mode) {mode_ = mode;}
  
  virtual int processCmd() throw(string);
};
//------------
//Process Calibration
//------------
class ProcessCalibration : public ProcessBase {
protected:
  string option_; //option of calibration procedure
  string  freqBAOCalibration_;//string MHz
  r_4 valfreqBAOCalibration_; //value MHz
  string bandWidthBAOCalibration_;//string MHz
  r_4 valbandWidthBAOCalibration_;//value MHz
  
  sa_size_t lowerFreqBin_;
  sa_size_t upperFreqBin_;

public:
  ProcessCalibration() {}
  virtual ~ProcessCalibration(){}

  void SetOption(const string& option) {option_ = option;}
  void SetFreqBAOCalibration(const string& freqBAOCalibration) { 
    freqBAOCalibration_ = freqBAOCalibration; 
    valfreqBAOCalibration_ = atof(freqBAOCalibration_.c_str());
  }
  void SetBandWidthBAOCalibration(const string& bandWidthBAOCalibration) { 
    bandWidthBAOCalibration_ = bandWidthBAOCalibration; 
    valbandWidthBAOCalibration_ = atof(bandWidthBAOCalibration_.c_str());
  }

  void ComputeLowerUpperFreqBin();
      
  virtual int processCmd() throw(string);
};
void ProcessCalibration::ComputeLowerUpperFreqBin() {
  sa_size_t c0 = round_sa(NUMBER_OF_FREQ*(valfreqBAOCalibration_-LOWER_FREQUENCY)/TOTAL_BANDWIDTH);
  sa_size_t dc = round_sa(NUMBER_OF_FREQ*valbandWidthBAOCalibration_/TOTAL_BANDWIDTH);
  lowerFreqBin_ = c0-dc/2;
  upperFreqBin_ = c0+dc/2;
}


//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{

  //Init process types
  map<string,IProcess*> process;
  process["dataOnOff"] = new ProcessONOFFData();
  process["gain"]      = new ProcessGain();
  process["calib"]     = new ProcessCalibration();

  //Init Sophya related modules
  //  SophyaInit();
  TArrayInitiator _inia; //nneded for TArray persistancy
  FitsIOServerInit(); //needed for input file

  //message used in Exceptions
  string msg;

  //Return code
  int rc = 0;

  //Arguments managements
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<11) return Usage(true);

  string action;
  string inputPath = "."; 
  string outputPath = "."; 
  string sourceName;
  string scaFile;
  string dateOfRun;
  string spectraDirectory;
  string freqBAOCalib = "";
  string bandWidthBAOCalib = "";
  string debuglev = "0";
  string mode = "";
  string numcycle;
  string calibrationOpt = "";  

  string typeOfFile="medfiltmtx";
  

  //  bool okarg=false;
  int ka=1;
  while (ka<(narg-1)) {
    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],"-source")==0) {
      sourceName=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-sca")==0) {
      scaFile=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-date")==0) {
      dateOfRun=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-specdir")==0) {
      spectraDirectory=arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-specname")==0) {
      typeOfFile=arg[ka+1];
      ka+=2;
    }    
    else if (strcmp(arg[ka],"-freqBAOCalib")==0) {
      freqBAOCalib = arg[ka+1];
      ka+=2;
    }
    else if (strcmp(arg[ka],"-bwBAOCalib")==0) {
      bandWidthBAOCalib = arg[ka+1];
      ka+=2;
    } 
    else if (strcmp(arg[ka],"-mode")==0) {
      mode =arg[ka+1];
      ka+=2; 
    }
    else if (strcmp(arg[ka],"-numcycle")==0) {
      numcycle =arg[ka+1];
      ka+=2; 
    }
    else if (strcmp(arg[ka],"-calibopt")==0) {
      calibrationOpt =arg[ka+1];
      ka+=2; 
    }
    else ka++;
  }//eo while


  //JEC 21/9/11 Give the input parameters START
  cout << " action " << action << "\n"
       << " inputPath = " << inputPath << "\n" 
       << " outputPath = " << outputPath << "\n"
       << " sourceName = " << sourceName << "\n"
       << " scaFile = " << scaFile << "\n"
       << " dateOfRun = " << dateOfRun << "\n"
       << " spectraDirectory = " << spectraDirectory << "\n"
       << " freqBAOCalib = " << freqBAOCalib  << "\n"
       << " bandWidthBAOCalib = " << bandWidthBAOCalib << "\n"
       << " debuglev = "  << debuglev  << "\n"
       << " mode = " << mode << "\n"
       << " numcycle = " << numcycle << "\n"
       << " calibrationOpt = " << calibrationOpt << endl;
  //JEC 21/9/11 Give the input parameters END


  try {
    //verification of action
    if(process.find(action) == process.end()) {
      msg = "action ";
      msg += action + " not valid... FATAL";
      rc = 999;
      throw msg;
    }
    




    //
    //Process initialisation...
    //
    try {
      ProcessBase* procbase = dynamic_cast<ProcessBase*>(process[action]);
      if (procbase == 0) {
        msg= "action ";
        msg += action + "Not a <process base> type...FATAL";
        rc = 999;
        throw msg;
      }
      procbase->SetInputPath(inputPath);
      procbase->SetOutputPath(outputPath);

      if ("" == sourceName) {
        msg = "(FATAL) missingsourceName  for action " + action;
        Usage(true);
        throw msg;
      }
      procbase->SetSourceName(sourceName);

      if ("" == dateOfRun) {
        msg = "(FATAL) missing dateOfRun for action " + action;
        Usage(true);
        throw msg;
      }
      procbase->SetDateOfRun(dateOfRun);

      
      if ("" == spectraDirectory) {
        msg = "(FATAL) missing spectraDirectory for action " + action;
        Usage(true);
        throw msg;
      }
      procbase->SetSpectraDirectory(spectraDirectory);

      if ("" == scaFile) {
        msg = "(FATAL) missing scaFile for action " + action;
        Usage(true);
        throw msg;
      }
      procbase->SetScaFileName(scaFile);

      if ("" == numcycle) {
        msg = "(FATAL) missing cycle number for action " + action;
        Usage(true);
        throw msg;
      }
      procbase->SetNumCycle(numcycle);


      procbase->SetTypeOfFile(typeOfFile);

      procbase->SetDebugLevel(debuglev);
    }
    catch(exception& e){
      throw e.what();
    }

    try {
      ProcessONOFFData* procdata = dynamic_cast<ProcessONOFFData*>(process[action]);
      if (procdata) {
        if (freqBAOCalib == "") {
          msg = "(FATAL) missing calibration BAO frequency for action " + action;
          Usage(true);
          throw msg;
        }
        procdata->SetFreqBAOCalibration(freqBAOCalib);
      }
    }
    catch(exception& e){
      throw e.what();
    }
    

    try {
      ProcessGain* procgain = dynamic_cast<ProcessGain*>(process[action]);
      if(procgain) {
        if (mode == "") {
          msg = "(FATAL) missing mode-type for action " + action;
          Usage(true);
          throw msg;
        }
        procgain->SetMode(mode);
      }
    }
    catch(exception& e){
      throw e.what();
    }

    try {
      ProcessCalibration* proccalib = dynamic_cast<ProcessCalibration*>(process[action]);
      if(proccalib) {
        if (calibrationOpt == "") {
          msg = "(FATAL) missing calibration option";
          Usage(true);
          throw msg;
        }
        if (freqBAOCalib == "") {
          msg = "(FATAL) missing calibration BAO frequency for action " + action;
          Usage(true);
          throw msg;
        }
        if (bandWidthBAOCalib == "") {
          msg = "(FATAL) missing calibration BAO frequency band width for action " + action;
          Usage(true);
          throw msg;
        }
        proccalib->SetOption(calibrationOpt);
        proccalib->SetFreqBAOCalibration(freqBAOCalib);
        proccalib->SetBandWidthBAOCalibration(bandWidthBAOCalib);
        proccalib->ComputeLowerUpperFreqBin();
      }
    }
    catch(exception& e){
      throw e.what();
    }

    //
    //execute command
    //
    rc = process[action]->processCmd();

  }
  catch (std::exception& sex) {
    cerr << "\n analyse.cc std::exception :"  << (string)typeid(sex).name() 
         << "\n msg= " << sex.what() << endl;
    rc = 78;
  }
  catch ( string str ) {
    cerr << "analyse.cc Exception raised: " << str << endl;
  }
  catch (...) {
    cerr << " analyse.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  


  cout << ">>>> analyse.cc ------- END ----------- RC=" << rc << endl;
  
  //Delete processes
  for_each(process.begin(),process.end(), DeleteMapFntor<string,IProcess*>());

  return rc;
}

//---------------------------------------------------
int Usage(bool flag) {
  cout << "Analyse.cc usage...." << endl;
  cout << "analyse  -act <action_type>: dataOnOff, gain, calib\n"
       << "         -inPath <path for input files: default='.'>\n"
       << "         -outPath <path for output files: default='.'>\n"
       << "         -source <source name> \n" 
       << "         -date <YYYYMMDD>\n"
       << "         -sca <file name scaXYZ.sum.trans>\n"
       << "         -specdir <generic directory name of spectra fits files>\n"
       << "         -specname <generic name of spectra fits files>\n"
       << "         -freqBAOCalib <freq in MHZ> freq. of calibration BAO\n"
       << "            valid for act=dataOnOff\n"
       << "         -bwBAOCalib <band width MHz> band width arround central freq. for calibration BAO\n"
       << "            valid for act=calib\n"
       << "         -mode <mode_type>:\n"
       << "            valid for act=gain, mode_type: On, Off\n"
       << "         -numcycle <number>,<number>:\n"
       << "            valid for all actions"
       << "         -calibopt <option>:\n"
       << "            valid for act=calib: indiv OR mean (NOT USED)"
       << "         -debuglev <number> [0=default]\n"
       << "           1: normal print\n"
       << "           2: save intermediate spectra\n"
       << endl;
  if (flag) {
    cout << "use <path>/analyse -h for detailed instructions" << endl;
    return 5;
  }
  return 1;
}

int ProcessBase::processCmd() throw(string) {
  int rc =0;
  string msg;
  if(debuglev_>0)cout << "Process Base" << endl;
  //------------------------
  //Use the sca file informations
  //------------------------
  //  string scaFullPathName = "./";
  //TOBE FIXED  scaFullPathName += sourceName_+"/"+dateOfRun_ + StringToLower(sourceName_)+"/";
  string scaFullPathName = inputPath_ + "/" 
    + sourceName_+ "/" +dateOfRun_ + StringToLower(sourceName_)+ "/" + scaFileName_;
  char* scaTupleColumnName[9] = {"cycle","stcalOn","spcalOn","stOn","spOn","stcalOff","spcalOff","stOff","spOff"};
  scaTuple_ = new NTuple(9,scaTupleColumnName);
  int n = scaTuple_->FillFromASCIIFile(scaFullPathName);
  if(n<0){ //Error
    msg = "(FATAL) NTuple error loading "+ scaFullPathName;
    rc = 999;
    throw msg;
  }
  
  if(debuglev_>1){
    cout << "ProcessBase::processCmd: dump tuple in " << scaFullPathName << endl;
    scaTuple_->Show(cout);
  }
  
  
  //Get the cycles (here consider consecutive cycles)    
  //The SCA file cannot be used as the DAQ can miss some cycles...
  //     r_8 firstCycle, lastCycle;
  //     scaTuple_->GetMinMax("cycle",firstCycle,lastCycle);
  //     ifirstCycle_ = (sa_size_t)firstCycle;
  //     ilastCycle_  = (sa_size_t)lastCycle;
  //Analyse the string given by -numcycle command line
  int ai1=0,ai2=0;
  sscanf(numcycle_.c_str(),"%d,%d",&ai1,&ai2);
  ifirstCycle_ = (sa_size_t)ai1;
  ilastCycle_  = (sa_size_t)ai2;
  

  //associate cycle number to index line in tuple
  sa_size_t nLines = scaTuple_->NbLines();
  for(sa_size_t iL=0; iL<nLines; ++iL){
    idCycleInTuple_[(sa_size_t)scaTuple_->GetCell(iL,"cycle")]=iL;
  }


  return rc;
}
//----------------------------------------------
int ProcessONOFFData::processCmd() throw(string) {
  int rc = 0;
  try {
    rc = ProcessBase::processCmd();
  } 
  catch (string s) {
    throw s;
  }
   if(debuglev_>0)cout << "Process Data" << endl;
  vector<string> modeList;
  modeList.push_back("On");
  modeList.push_back("Off");
  vector<string>::const_iterator iMode;
  
  uint_4 id; 
  string tag;

  //
  //Process to get sucessively
  //Raw Spectra, 
  //BAO Calibrated Spectra 
  //and RT Calibrated Spectra
  //The pocesses are separated to allow intermediate save of results

  map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
  map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
  
  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
    string mode = *iMode;
    if(debuglev_>0)cout << "Process RAW Mode " << mode << endl;

    //------------------------------------------
    //Produce Raw spectra per cycle
    //------------------------------------------

    string directoryName;
    list<string> listOfSpecFiles;
    list<string>::const_iterator iFile, iFileEnd;
    
        
    //
    //loop on cycles
    //
    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
      //TOBE FIXED      directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
      directoryName = "./" + mode + "/";
      stringstream sicycle;
      sicycle << icycle;
      directoryName += spectraDirectory_ + sicycle.str() + "/";

      //read directory
      listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
      

      //compute mean of spectra created in a cycle
      if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
      TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
      iFileEnd = listOfSpecFiles.end();
      r_4 nSpectres  = 0;
      for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
        aSpectrum >> spectre;
        spectreMean += spectre;
        nSpectres++;
      }// end of for files 
      if(nSpectres>0) spectreMean /= nSpectres;
      
      //save mean spectrum
      spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
    }//end of for cycles
  }//end loop on mode for raw preocess

  if(debuglev_>1) {//save mean spectra on file 
    cout << "Save mean raw spectra" << endl;
    string fileName;
    //TOBE FIXED    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataRaw" + ".ppf";
    fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";

    POutPersist fos(fileName);
    id=0;
    iSpectreEnd = spectreCollect.end();
    for (iSpectre = spectreCollect.begin();
         iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
      tag = "specRaw";

      //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
//       stringstream sid;
//       sid << id;
//       tag += sid.str();
      tag += (iSpectre->first).first;
      stringstream sid;
      sid << (iSpectre->first).second;
      tag += sid.str();
      if(debuglev_>9) {
        cout << "save tag<" << tag << ">" << endl;
      }
      //JEC 20/9/11 modif tag to take into account Mode and Cycle number END

      fos << PPFNameTag(tag) << iSpectre->second;
    }
  }//end of save fits
  


  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
    string mode = *iMode;
    if(debuglev_>0)cout << "Process CALIB BAO Mode " << mode << endl;
    //------------------------------------------
    // Correct Raw spectra for BAO calibration
    //------------------------------------------
    //Read BAO calibration files
    sa_size_t nr,nc; //values read 
//JEC 20/9/11 use mean calibration coeff upon all cycles START
//    bool first = true;
//    TMatrix<r_4> calibBAOfactors;
//     for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
//       string calibFileName = inputPath_+ "/" 
//      + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) 
//      + "/calib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
//      + mode + "_" + freqBAOCalibration_ + "MHz";
//       stringstream channel;
//       channel << iCh;
//       calibFileName += "-Ch" + channel.str() + "Cycles.txt";
//       if(debuglev_>0) cout << "Read file " << calibFileName << endl;
//       ifstream ifs(calibFileName.c_str());
//       if ( ! ifs.is_open() ) {
//      rc = 999;
//      throw calibFileName + " cannot be opened...";
//       }      
//       TVector<r_4> aCalib;
//       if(debuglev_>9) cout << "Debug 1" << endl;
//       aCalib.ReadASCII(ifs,nr,nc);
//       if(debuglev_>9) cout << "Debug 2" << endl;
//       if(first){
//      first = false;
//      if(debuglev_>9) cout << "Debug 3" << endl;
//      calibBAOfactors.SetSize(NUMBER_OF_CHANNELS,nr);
//      if(debuglev_>9) cout << "Debug 4" << endl;
//       }
//       calibBAOfactors( Range(iCh), Range::all() ) = aCalib.Transpose();
//       if(debuglev_>9) cout << "Debug 5" << endl;
//       ifs.close();
//     }//end of loop on channels

    
    string calibFileName = inputPath_+ "/" 
      + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) 
      + "/calib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
      + mode + "_" + freqBAOCalibration_ + "MHz-All.txt";

    if(debuglev_>0) cout << "Read file " << calibFileName << endl;
    ifstream ifs(calibFileName.c_str());
    if ( ! ifs.is_open() ) {
      rc = 999;
      throw calibFileName + " cannot be opened...";
    }   
    TVector<r_4> calibBAOfactors;
    if(debuglev_>9) cout << "Debug 1" << endl;
    calibBAOfactors.ReadASCII(ifs,nr,nc);
    if(debuglev_>9){
      cout << "Debug 2: (nr,nc): "<< nr << "," << nc << endl;
      calibBAOfactors.Print(cout);
    }

//JEC 20/9/11 use mean calibration coeff upon all cycles END

    //
    //spectra corrected by BAO calibration factor 
    //-----make it different on Channels and Cycles (1/06/2011) OBSOLETE
    //use mean calibration coeff upon all cycles (20/6/11)
    //warning cycles are numbered from 1,...,N
    //
    if(debuglev_>0)cout << "do calibration..." << endl;
    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
      TMatrix<r_4> specmtx(spectreCollect[make_pair(mode,icycle)],true); //share the memory 
      
      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
        //JEC 20/9/11 use mean calibration coeff upon all cycles START

        //      specmtx( Range(iCh), Range::all() ) /= calibBAOfactors(iCh,icycle-1);
        specmtx( Range(iCh), Range::all() ) /= calibBAOfactors(iCh);
        //JEC 20/9/11 use mean calibration coeff upon all cycles END
      }
    }
  } //end loop mode for BAO calib

  if(debuglev_>1){ //save mean spectra BAO calibrated on file
    cout << "save calibrated BAO spectra" << endl;
    string fileName;
    //TO BE FIXED    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAOCalib" + ".ppf";
    fileName = "./dataBAOCalib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";

    POutPersist fos(fileName);
    iSpectreEnd = spectreCollect.end();
    id=0;
    for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {

      tag = "specBAOCalib";
      //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
//       stringstream sid;
//       sid << id;
//       tag += sid.str();
      tag += (iSpectre->first).first;
      stringstream sid;
      sid << (iSpectre->first).second;
      tag += sid.str();
      if(debuglev_>9) {
        cout << "save tag<" << tag << ">" << endl;
      }
      //JEC 20/9/11 modif tag to take into account Mode and Cycle number END

      fos << PPFNameTag(tag) << iSpectre->second;
    }
  }//end of save fits

  
  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
    string mode = *iMode;
    if(debuglev_>0)cout << "Process CALIB RT Mode " << mode << endl;
    //------------------------------------------
    // Correct BAO calib spectra for RT calibration
    //------------------------------------------
    //Very Preliminary May-June 11
    //coef RT @ 1346MHz Ouest - Est associees a Ch 0 et 1
    r_4 calibRT[NUMBER_OF_CHANNELS] = {27.724, 22.543};
    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
      TMatrix<r_4> specmtx(spectreCollect[make_pair(mode,icycle)],true); //share the memory    
      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
        specmtx( Range(iCh), Range::all() ) *= calibRT[iCh];
      }
    }
  }//end loop on mode RT calib

  {//save mean spectra BAO & RT calibrated on file
    if(debuglev_>0)cout << "save calibrated BAO & RT spectra" << endl;
    string fileName;
    //TO BE FIXED    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAORTCalib" + ".ppf";
     fileName = "./dataBAORTCalib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";   
    POutPersist fos(fileName);
    iSpectreEnd = spectreCollect.end();
    id = 0;
    for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
      tag = "specBAORTCalib";
      //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
//       stringstream sid;
//       sid << id;
//       tag += sid.str();
      tag += (iSpectre->first).first;
      stringstream sid;
      sid << (iSpectre->first).second;
      tag += sid.str();
      if(debuglev_>9) {
        cout << "save tag<" << tag << ">" << endl;
      }
      //JEC 20/9/11 modif tag to take into account Mode and Cycle number END
      fos << PPFNameTag(tag) << iSpectre->second;
    }
  }//end of save fits

  //------------------------------------------
  // Perform ON-OFF
  //------------------------------------------
  
  map< sa_size_t, TMatrix<r_4> > diffCollect;
  map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd;

  TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero
  r_4 nCycles = 0;
  for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
    nCycles++;
    TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory 
    TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory 
    TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff;
    diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) ));
    diffMeanOnOff += diffOnOff;
  }//end loop on cycle
  if(nCycles>0) diffMeanOnOff/=nCycles;

  //exctract channels and do the mean
  TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0
  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
    meanOfChan += diffMeanOnOff.Row(iCh).Transpose();
  }
  meanOfChan /= (r_4)NUMBER_OF_CHANNELS;
  


  {//save diff ON-OFF BAO & RT calibrated
    if(debuglev_>0)cout << "save ON-OFF spectra" << endl;
    string fileName;
    //TO BE FIXED    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "diffOnOff" + ".ppf";
    fileName = "./diffOnOff_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
    POutPersist fos(fileName);
    iDiffEnd = diffCollect.end();
    id = 0;
    for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
      tag = "specONOFF";
      stringstream sid;
      //JEC 20/9/11      sid << id;
      sid << iDiff->first;
      tag += sid.str();
      if(debuglev_>9) {
        cout << "save tag<" << tag << ">" << endl;
      }
      fos << PPFNameTag(tag) << iDiff->second;
    }
    //save the mean also
    fos << PPFNameTag("specONOFFMean") << diffMeanOnOff;
    fos << PPFNameTag("specONOFF2ChanMean") << meanOfChan;
  }//end of save fits

  return rc;
} //ProcessONOFFData::processCmd
//
//----------------------------------------------
//
int ProcessGain::processCmd() throw(string) {
  int rc = 0;
  string msg = "";

  try {
    rc = ProcessBase::processCmd();
  } 
  catch (string s) {
    throw s;
  }
  if(debuglev_>0)cout << "Process Gain" << endl;
  
  string directoryName;
  //TOBE FIXED  directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode_ + "/";
  //JEC 6/09/2011 numcycle_ repalced by ifirstCycle_ according to definition of numcycle_ and the fact that for GAIN 1 cycle is involved
  stringstream thegaincycle;
  thegaincycle << ifirstCycle_;
  directoryName = inputPath_ + "/" 
    + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/" 
    + mode_ + "/" + spectraDirectory_ + thegaincycle.str()  + "/";
  
  list<string> listOfSpecFiles;
  list<string>::const_iterator iFile, iFileEnd;
  //read directory

  listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
  
  //Mean power computed over the N channels to select the spectra for gain computation 
  //The threshold is computed "on-line" as  1%  of the difference (max power -min power) over the
  // the min power.
  //A possible alternative is to set an absolute value...
  if(debuglev_>0)cout << "compute mean poser spectra for files in " << directoryName << endl;
  iFileEnd = listOfSpecFiles.end();
  map<string, r_4> meanpowerCollect;
  //JEC 21/9/11 add meanpower for each Channels START
  map<string, r_4> meanPowerPerChanCollect;
  //JEC 21/9/11 add meanpower for each Channels END

  map<string, r_4>::const_iterator iMeanPow, iMeanPowEnd;

  for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
    FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
    TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
    aSpectrum >> spectre;
    meanpowerCollect[*iFile] = spectre.Sum()/spectre.Size();
    //JEC 21/9/11 add meanpower for each Channels START
    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
      TVector<r_4> specChan(NUMBER_OF_FREQ);
      specChan = spectre.Row(iCh).Transpose();
      stringstream tmp;
      tmp << iCh;
      string tag = *iFile + "_" + tmp.str();
      meanPowerPerChanCollect[tag] = specChan.Sum()/specChan.Size();
    }
    //JEC 21/9/11 add meanpower for each Channels END
  }//end of for files
  pair<string, r_4> minelement = *min_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
  pair<string, r_4> maxelement = *max_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
  if(debuglev_>1){
    cout << "meanpowerCollect has " << meanpowerCollect.size() << " spectra registered" << endl;
    cout << "find min mean power "<<minelement.second << " for " << minelement.first << endl;
    cout << "find max mean power "<<maxelement.second << " for " << maxelement.first << endl;
  }
  r_4 threshold = minelement.second + 0.01*(maxelement.second-minelement.second);
  if(debuglev_>1){
    cout << "threshold found at " << threshold <<endl;
  }  

  TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
  r_4 nSpectres  = 0;
  iMeanPowEnd = meanpowerCollect.end();
  for (iMeanPow = meanpowerCollect.begin(); iMeanPow != iMeanPowEnd; ++iMeanPow) {
    if ( iMeanPow->second <= threshold ) {
      //TODO avoid the reloading of the file may speed up
      FitsInOutFile aSpectrum(iMeanPow->first,FitsInOutFile::Fits_RO);
      TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
      aSpectrum >> spectre;
      spectreMean += spectre;
      nSpectres++;
    }
  }
  if(debuglev_>1){
    cout << "Number of files selected for gain " << nSpectres <<endl;
  }    
  if(nSpectres>0) {
    spectreMean /= nSpectres;
  } else {
    stringstream tmp;
    tmp << threshold;
    msg = "Gain: cannot find a spectra above threshold value =" + tmp.str() + " ... FATAL";
    throw msg;
  }

  //Save gain spectra
  {
    //use ! to override the file: special features of cfitsio library
    string fileName;
    //TOBE FIXED   fileName = "!./" + sourceName_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
    fileName = "!"+ outputPath_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
    if(debuglev_>1){
      cout << "save gains in " << fileName << endl;
    }
    FitsInOutFile fos(fileName, FitsInOutFile::Fits_Create);
    fos << spectreMean;
  }
  //save mean power values
  {
    vector<r_4> tmp;
    //JEC 21/9/11 add meanpower for each Channels START
    vector<r_4> tmpCh0; //for Chan 0
    vector<r_4> tmpCh1; //for Chan 1
    //JEC 21/9/11 add meanpower for each Channels END
    for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
      if (debuglev_>9) {
        cout << "Gain: save mean power of file: " << *iFile << endl;
      }
      tmp.push_back(meanpowerCollect[*iFile]);
      //JEC 21/9/11 add meanpower for each Channels START  
      stringstream tmp0;
      tmp0 << (sa_size_t)0;
      string tag0 = *iFile + "_" + tmp0.str();
      tmpCh0.push_back(meanPowerPerChanCollect[tag0]);
      if (NUMBER_OF_CHANNELS>1){
        stringstream tmp1;
        tmp1 << (sa_size_t)1;
        string tag1 = *iFile + "_" + tmp1.str();
        tmpCh1.push_back(meanPowerPerChanCollect[tag1]);
      }
      //JEC 21/9/11 add meanpower for each Channels END
    }
    string fileName;
    //TOBE FIXED    fileName = "./" + sourceName_ + "/gain_monitor_" + dateOfRun_ 
    fileName = outputPath_ + "/gain_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
    POutPersist fos(fileName); 
    TVector<r_4> tvtmp(tmp);
    fos.PutObject(tvtmp,"gainmoni"); //Attention initialement le tag etait "monitor"...
    //JEC 21/9/11 add meanpower for each Channels START  
    TVector<r_4> tvtmp0(tmpCh0);
    fos.PutObject(tvtmp0,"gainmoni0");
    if (NUMBER_OF_CHANNELS>1){
      TVector<r_4> tvtmp1(tmpCh1);
      fos.PutObject(tvtmp1,"gainmoni1");
    }
    //JEC 21/9/11 add meanpower for each Channels END
  }
  
  cout << "OK gain finished" <<endl;
   return rc;
}//ProcessGain::processCmd
//
//----------------------------------------------
//
int ProcessCalibration::processCmd() throw(string) {
  int rc = 0;
  string msg = "";

  try {
    rc = ProcessBase::processCmd();
  } 
  catch (string s) {
    throw s;
  }
  if(debuglev_>0)cout << "Process Calibration with option "<< option_ << endl;
  
  vector<string> modeList;
  modeList.push_back("On");
  modeList.push_back("Off");

  r_8 t0absolute = -1;  //TIMEWIN of the first spectra used
  r_8 timeTag0   = -1;   //MEANTT, mean TIME TAG of the first paquet window  
  bool first = true;     // for initialisation
  
  // Power Tuple
  // mode, cycle, time, {Ch0, ... ,ChN} mean poser in the range [f0-Bw/2, f0+Bw/2] 
  vector<string> varPowerTupleName; //ntuple variable naming
  varPowerTupleName.push_back("mode");
  varPowerTupleName.push_back("cycle");  
  varPowerTupleName.push_back("time");
  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
    stringstream tmp;
    tmp << iCh;
    varPowerTupleName.push_back("Ch"+tmp.str());
  }
  r_4 valPowerTuple[varPowerTupleName.size()];
  NTuple powerEvolution(varPowerTupleName); 
  
  
  //-----------------
  //Start real process
  //------------------
  for (size_t iMode = 0; iMode < modeList.size(); ++iMode) {
    string mode = modeList[iMode];
    if(debuglev_>0)cout << "Process Calibration for Mode " << mode << endl;

    //initialize the start of each calibration procedure given by the RT SCA file
    //see ProcessBase::processCmd for the structure of the sca file
    string scaStartCalibName = "stcal"+mode; 
    sa_size_t idStartCalib   = scaTuple_->ColumnIndex(scaStartCalibName);

    string directoryName;
    list<string> listOfSpecFiles;
    list<string>::const_iterator iFile, iFileEnd;            
    //
    //loop on cycles
    //
    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {

      directoryName = inputPath_ + "/" 
        + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
      stringstream sicycle;
      sicycle << icycle;
      directoryName += spectraDirectory_ + sicycle.str() + "/";
      
      //read directory
      listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);

      iFileEnd = listOfSpecFiles.end();
      DVList header;
      TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0

      for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
        aSpectrum.GetHeaderRecords(header);
        aSpectrum >> spectre;

        if(first){//initialise the timer
          first = false;
          t0absolute = header.GetD("TIMEWIN")/1000.;
          timeTag0 = header.GetD("MEANTT");
          if (debuglev_>1){
            cout << "debug Header of " << *iFile << endl;
            cout << "TIMEWIN = " << setprecision(12) << t0absolute << " "
                 << "MEANTT = " << setprecision(12) << timeTag0 << endl;
          }
        }//end init timer
        
        //Set time of current file
        //Add to the non-precise absolute origin an precise increment
        r_4 timeTag = t0absolute + (header.GetD("MEANTT") - timeTag0);
        int index=0;
        valPowerTuple[index++] = iMode;
        valPowerTuple[index++] = icycle;
        valPowerTuple[index++] = timeTag-scaTuple_->GetCell(idCycleInTuple_[icycle],idStartCalib); //take the RT time start as refernce 

        //Compute the mean power of the two channel (separatly) in the range [f-bw/2, f+bw/2]
        for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
          TMatrix<r_4> tmp(spectre(Range(iCh),Range(lowerFreqBin_,upperFreqBin_)),false);
          r_4 mean = tmp.Sum()/(r_4)tmp.Size();
          valPowerTuple[index++] = mean;
        }//end of channel loop
       
        //Fill Tuple
        powerEvolution.Fill(valPowerTuple);
      }//end of files loop
    }//end of cycles loop
  }//end of mode loop

  //store power evolution Tuple
  if(debuglev_>0){
    cout << "ProcessCalibration::processCmd: dump powerEvolution tuple" << endl;
    powerEvolution.Show(cout);
  }
  //
  //Compute Calibration Coefficients as function of mode/cycle/channels
  //

  //Tsys,Incremant Intensity... Tuple
  //mode mode cycle [(tsys0,dint0),...,(tsysN,dintN)]
  vector<string> varCalibTupleName;
  varCalibTupleName.push_back("mode");
  varCalibTupleName.push_back("cycle");
  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
    stringstream tmp;
    tmp << iCh;
    varCalibTupleName.push_back("tsys"+tmp.str());
    varCalibTupleName.push_back("dint"+tmp.str());
  }
  r_4 valCalibTuple[varCalibTupleName.size()];
  NTuple calibEvolution(varCalibTupleName);

  // select time E [twin0,twin1] U [twin2,twin3] for Tsys
  // time unit = second 
  const r_4 twin0 = -3.;
  const r_4 twin1 = -1.;
  const r_4 twin2 =  6.;
  const r_4 twin3 =  8;
  const r_4 thresholdFactor = 0.80; //80% of the diff. Max-Min intensity

  sa_size_t inModeIdx = powerEvolution.ColumnIndex("mode");
  sa_size_t inCycleIdx= powerEvolution.ColumnIndex("cycle");
  sa_size_t inTimeIdx = powerEvolution.ColumnIndex("time");
  sa_size_t inOffsetCh0 = powerEvolution.ColumnIndex("Ch0"); //nb Ch0 position in the powerEvolution tuple  
  if(debuglev_>1) cout << "DEBUG: in Idx: (" 
                       << inModeIdx << ", "
                       << inCycleIdx << ", "
                       << inTimeIdx << ", "
                       << inOffsetCh0 << ")"
                       << endl;

 
  size_t outModeIdx = calibEvolution.ColumnIndex("mode");
  size_t outCycleIdx= calibEvolution.ColumnIndex("cycle");
  size_t outOffsetCh0 = calibEvolution.ColumnIndex("tsys0"); //nb Ch0 position in the monitor tuple  
  size_t outNDataPerCh= 2;
  if(debuglev_>1)  cout << "DEBUG: out Idx: (" 
                        << outModeIdx << ", "
                        << outCycleIdx << ", "
                        << outOffsetCh0 << ")"
                        << endl;

  sa_size_t nPowerEvolEntry = powerEvolution.NEntry();
  double minMode;
  double maxMode;
  powerEvolution.GetMinMax("mode",minMode,maxMode);
  double minCycleNum;
  double maxCycleNum;
  powerEvolution.GetMinMax("cycle",minCycleNum,maxCycleNum);
  if(debuglev_>1)  cout << "DEBUG mode ("<<minMode<<","<<maxMode<<")\n"
                        << "cycle ("<<minCycleNum<<","<<maxCycleNum<<")"
                        << endl;

  r_4* aPowerEvolEntry; // a ligne of the powerEvolution tuple
//   r_4* aPowerEvolEntry = new r_4[powerEvolution.NbColumns()]; // a ligne of the powerEvolution tuple

  r_4 minMean;
  r_4 maxMean;

  for (size_t iMode = (size_t)minMode; iMode <= (size_t)maxMode; iMode++){
    for (size_t iCycle = (size_t)minCycleNum; iCycle <= (size_t)maxCycleNum; iCycle++ ){
      if(debuglev_>1) cout<<"DEBUG >>>>>>> mode/cycle: " << iMode << "/" << iCycle << endl;
 
      valCalibTuple[outModeIdx]=iMode;
      valCalibTuple[outCycleIdx]=iCycle;
      //
      //Compute the Min && Max value of each Ch<i> data
      if(debuglev_>1) cout<<"DEBUG compute Min and Max for each channels" << endl;
      //
      TVector<r_4> chMean[NUMBER_OF_CHANNELS];
      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
        chMean[iCh].SetSize(nPowerEvolEntry);
      }
      for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
        aPowerEvolEntry = powerEvolution.GetVec(ie);
        if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
          for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
            chMean[iCh](ie) = aPowerEvolEntry[iCh+inOffsetCh0];
          }//end cut
        }
      }//eo loop on tuple entries
      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
        //
        //select values to compute background Tsys
        if(debuglev_>1) {
          cout<< "DEBUG >>>> Ch[" << iCh << "]" << endl;
          cout<< "DEBUG select values to compute background Tsys " << endl;
        }
        //
       
        std::vector<r_4> lowerInt;
        for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
          aPowerEvolEntry = powerEvolution.GetVec(ie);
          if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
            r_4 time = aPowerEvolEntry[inTimeIdx];
            if ( (time >= twin0 && time <= twin1) ||
                 (time >= twin2 && time <= twin3)
                 ) lowerInt.push_back(aPowerEvolEntry[iCh+inOffsetCh0]);
          }//end cut
        } //end loop entries
        //
        //compute the Tsys
        if(debuglev_>1) cout <<"DEBUG compute Tsys" << endl;
        //
        std::nth_element(lowerInt.begin(),
                         lowerInt.begin()+lowerInt.size()/2,
                         lowerInt.end());
        r_4 tsys = *(lowerInt.begin()+lowerInt.size()/2);
        if(debuglev_>1) cout << "DEBUG tsys["<<iCh<<"] : " << tsys <<endl;
        //
        //set the threshold for DAB detection
        //
        chMean[iCh].MinMax(minMean,maxMean);
        minMean = (tsys > minMean) ? tsys : minMean; //pb of empty vector
        if(debuglev_>1) cout << "DEBUG min["<<iCh<<"] : " << minMean
                             << " max["<<iCh<<"] : " << maxMean
                             <<endl;
        r_4 deltathres = thresholdFactor * (maxMean-minMean);
        r_4 thres =  tsys + deltathres;
        if(debuglev_>1) cout << "DEBUG thres["<<iCh<<"] : " << thres <<endl;
        //
        //collect upper part of the DAB
        if(debuglev_>1) cout <<"DEBUG collect upper part of the DAB" << endl;
        //
        std::vector<r_4> upperInt;
        for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
          aPowerEvolEntry = powerEvolution.GetVec(ie);
          if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
            r_4 mean = aPowerEvolEntry[iCh+inOffsetCh0];
            if (mean >= thres) upperInt.push_back(mean);
          }//end cut
        }//eo loop on channels
        //
        //compute adjacent differences to detect the 2 DAB levels
        //
        if(debuglev_>1) cout << "(DEBUG )size upper [" << iCh << "]: " << upperInt.size() <<endl;
        std::vector<r_4> upperIntAdjDiff(upperInt.size());
        std::adjacent_difference(upperInt.begin(),
                                 upperInt.end(),
                                 upperIntAdjDiff.begin());
        //
        //Search the 2 DAB states
        if(debuglev_>1) cout<<"DEBUG Search the 2 DAB states" << endl;
        //
        std::vector<r_4> upIntState[2];
        int state=-1;
        for (size_t i=1;i<upperInt.size();i++) {//skip first value
          if (fabs(upperIntAdjDiff[i])<upperInt[i]*0.010) {
            if(state == -1) state=0;
            if(state == -2) state=1;
            upIntState[state].push_back(upperInt[i]);
          } else {
            if (state == 0) state=-2;
          }
        }
        //
        //Take the mean of the median values of each levels
        if(debuglev_>1)cout << "DEBUG mean of the median values of each levels" << endl;       
        //
        r_4 meanUpper = 0;
        r_4 nval = 0;
        for (int i=0;i<2;i++) {
          if (!upIntState[i].empty()) {
            std::nth_element(upIntState[i].begin(),
                             upIntState[i].begin()+upIntState[i].size()/2,
                             upIntState[i].end());
            meanUpper += *(upIntState[i].begin()+upIntState[i].size()/2);
            nval++;
          } 
        }
        meanUpper /= nval;
        //
        //Finaly the increase of power due to the DAB is:
        //
        r_4 deltaInt = meanUpper - tsys;
        if(debuglev_>1) cout << "DEBUG deltaInt["<<iCh<<"] : " << deltaInt <<endl;
        //
        //Save for monitoring and final calibration operations
        //
        valCalibTuple[outOffsetCh0+outNDataPerCh*iCh]   = tsys;
        valCalibTuple[outOffsetCh0+outNDataPerCh*iCh+1] = deltaInt;
      }//end loop on channels
      if(debuglev_>1) cout<<"DEBUG Fill the calibEvolution tuple" << endl;
      calibEvolution.Fill(valCalibTuple);
    }//eo loop on Cycles
  }//eo loop on Modes
  //
  //store calibration evolution Tuple
  //
  if(debuglev_>0){
    cout << "ProcessCalibration::processCmd: dump calibEvolution tuple" << endl;
    calibEvolution.Show(cout);
  }
  //
  //Compute the means per channel of the calibration coefficients (deltaInt)
  //and save cycle based Calibration Ctes in file named as 
  //    <source>-<date>-<mode>-<fcalib>MHz-Ch<channel>Cycles.txt
  //   format <cycle> <coef> 
  //the means are stored in files     
  //    <source>-<date>-<mode>-<fcalib>MHz-All.txt
  //   format <channel> <coef>
  //
  sa_size_t nModes = (sa_size_t)(maxMode - minMode) + 1;
  string calibCoefFileName;
  ofstream  calibMeanCoefFile[nModes]; //Mean over cycles
  ofstream  calibCoefFile[nModes][NUMBER_OF_CHANNELS]; //cycle per cycle
  for (sa_size_t iMode=0; iMode<nModes; iMode++){
    //The file for the Means Coeff.
    calibCoefFileName = outputPath_ + "/calib_" 
      + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
      + modeList[iMode] + "_"
      + freqBAOCalibration_ + "MHz-All.txt";
    calibMeanCoefFile[iMode].open(calibCoefFileName.c_str());
    //The file for the cycle per cycle Coeff.
    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
      stringstream chString;
      chString << iCh;
      calibCoefFileName = outputPath_ + "/calib_" 
        + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
        + modeList[iMode] + "_"
        + freqBAOCalibration_ + "MHz-Ch" + chString.str() + "Cycles.txt";
      calibCoefFile[iMode][iCh].open(calibCoefFileName.c_str());
    }
  }

  r_4* aCalibEvolEntry;
  sa_size_t nCalibEvolEntry = calibEvolution.NEntry();

  TMatrix<r_4> meanCalibCoef(nModes,NUMBER_OF_CHANNELS); //by default init to 0
  TMatrix<r_4> nData4Mean(nModes,NUMBER_OF_CHANNELS);

  //READ the calibration tuple, fill the array for mean and write to ascii file
  for (sa_size_t ie=0; ie<nCalibEvolEntry; ie++){
    aCalibEvolEntry = calibEvolution.GetVec(ie);
    if(debuglev_>1){
      cout << "DEBUG mode/cycle/Coef " 
           << aCalibEvolEntry[outModeIdx] << " "
           << aCalibEvolEntry[outCycleIdx] << " ";
    }
    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
      if(debuglev_>1) cout << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1] << " "; // for debug
      sa_size_t currentMode = (sa_size_t)(aCalibEvolEntry[outModeIdx] - minMode); //ok even if minMode <> 0
      nData4Mean(currentMode,iCh)++;
      meanCalibCoef(currentMode,iCh) += aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1];

      calibCoefFile[currentMode][iCh] << aCalibEvolEntry[outCycleIdx] << " "
                                      << setprecision(12)
                                      << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1]
                                      << endl;
    }
    if(debuglev_>1) cout << endl; //for debug
  }
  
  //Perform means: true means that div per 0 is treated (slower but safer) 
  meanCalibCoef.Div(nData4Mean,true);
  
  if(debuglev_>0){
    cout << "DEBUG nData4Mean" << endl;
    nData4Mean.Print(cout);
    cout << "DEBUG meanCalibCoef (averaged)" << endl;
    meanCalibCoef.Print(cout);
  }

  //Save means Coef
  for (sa_size_t iMode=0; iMode<nModes; iMode++){
    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
      calibMeanCoefFile[iMode] << setprecision(12)
                               << meanCalibCoef(iMode,iCh)
                               << endl;
    }
  }

  //Save Monitor File
  {
    string fileName = outputPath_ + "/calib_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
    POutPersist calibMonitorFile(fileName);
    calibMonitorFile << PPFNameTag("powermoni") <<  powerEvolution;
    calibMonitorFile << PPFNameTag("calibmoni") <<  calibEvolution;
  }

  //Clean & Return
  for (sa_size_t iMode=0; iMode<nModes; iMode++){
    calibMeanCoefFile[iMode].close();
    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
      calibCoefFile[iMode][iCh].close();
    }
  }

  cout << "Ok calibration finished" << endl;  
  return rc;
}//ProcessCalibration::processCmd