//
// 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 <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 uint_4 NUMBER_OF_CHANNELS = 2;
const uint_4 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 (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
}
//-----
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;
  }
  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 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 SetScaFileName(const string& scaFileName) { scaFileName_ =scaFileName; }

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

  virtual int processCmd() throw(string);
  
protected:
  string sourceName_;
  string dateOfRun_;
  string spectraDirectory_;
  string typeOfFile_;
  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
  string numcycle_; //cycle number
  uint_4 inumcycle_;// "
public:
  ProcessGain(){}
  virtual ~ProcessGain(){}

  void SetMode(const string& mode) {mode_ = mode;}
  void SetNumCycle(const string& numcycle) {
    numcycle_ = numcycle;
    inumcycle_ = atoi(numcycle_.c_str());
  }
  
  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 sourceName;
  string scaNumber;
  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;
    }
    if (strcmp(arg[ka],"-act")==0) {
      action=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) {
      scaNumber=arg[ka+1];
      scaFile = "sca"+scaNumber;
      scaFile += ".sum.trans";
      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



  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->SetSourceName(sourceName);
      procbase->SetDateOfRun(dateOfRun);
      procbase->SetSpectraDirectory(spectraDirectory);
      procbase->SetScaFileName(scaFile);
//       procbase->SetFirstCycle((uint_4)firstCycle);
//       procbase->SetLastCycle((uint_4)lastCycle);
      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);
	if (numcycle == "") {
	  msg = "(FATAL) missing cycle number for action " + action;
	  Usage(true);
	  throw msg;
	}
	procgain->SetNumCycle(numcycle);
      }
    }
    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>: data, gain, calib\n"
       << "         -source <source name> \n" 
       << "         -date <YYYYMMDD>\n"
       << "         -sca <number XYZ in the scaXYZ file>\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=data\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"
       << "         -calibopt <option>:\n"
       << "            valid for act=calib: indiv OR mean"
       << "         -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 = "./";
  scaFullPathName += sourceName_+"/";
  scaFullPathName += 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;
  }
  
  scaTuple_->Show(cout); //debug
  
  //Get the cycles (here consider consecutive cycles)
  r_8 firstCycle, lastCycle;
  scaTuple_->GetMinMax("cycle",firstCycle,lastCycle);
  ifirstCycle_ = (uint_4)firstCycle;
  ilastCycle_  = (uint_4)lastCycle;

  //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++) {
      directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +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;
    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataRaw" + ".ppf";
    POutPersist fos(fileName);
    id=0;
    iSpectreEnd = spectreCollect.end();
    for (iSpectre = spectreCollect.begin();
	 iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
      tag = "specRaw";
      stringstream sid;
      sid << id;
      tag += sid.str();
      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
    TMatrix<r_4> calibBAOfactors;
    sa_size_t nr,nc; //values read 
    bool first = true;

    for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
      string calibFileName = "./" + sourceName_ + "/" 
	+ 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;
      aCalib.ReadASCII(ifs,nr,nc);
      if(first){
	first = false;
	calibBAOfactors.SetSize(NUMBER_OF_CHANNELS,nr);
      }
      calibBAOfactors( Range(iCh), Range::all() ) = aCalib.Transpose();
      ifs.close();
    }//end of loop on channels


    //
    //spectra corrected by BAO calibration factor 
    //make it different on Channels and Cycles (1/06/2011)
    //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){
	specmtx( Range(iCh), Range::all() ) /= calibBAOfactors(iCh,icycle-1);
      }
    }
  } //end loop mode for BAO calib

  if(debuglev_>1){ //save mean spectra BAO calibrated on file
    cout << "save calibrated BAO spectra" << endl;
    string fileName;
    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAOCalib" + ".ppf";
    POutPersist fos(fileName);
    iSpectreEnd = spectreCollect.end();
    id=0;
    for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
      tag = "specBAOCalib";
      stringstream sid;
      sid << id;
      tag += sid.str();
      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;
    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAORTCalib" + ".ppf";
    POutPersist fos(fileName);
    iSpectreEnd = spectreCollect.end();
    id = 0;
    for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
      tag = "specBAORTCalib";
      stringstream sid;
      sid << id;
      tag += sid.str();
      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;
    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "diffOnOff" + ".ppf";
    POutPersist fos(fileName);
    iDiffEnd = diffCollect.end();
    id = 0;
    for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
      tag = "specONOFF";
      stringstream sid;
      sid << id;
      tag += sid.str();
      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;
  directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode_ + "/";
  directoryName += spectraDirectory_ + numcycle_ + "/";
  
  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;
  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();
  }//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;
    fileName = "!./" + sourceName_ + "/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;
    for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
      tmp.push_back(meanpowerCollect[*iFile]);
    }
    string fileName;
    fileName = "./" + sourceName_ + "/gain_monitor_" + dateOfRun_ 
      + "_" + StringToLower(sourceName_) + ".ppf";
    POutPersist fos(fileName); 
    TVector<r_4> tvtmp(tmp);
    fos.PutObject(tvtmp,"monitor");
  }

   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
  
  //Tuple
  // mode, cycle, time, {Ch0, ... ,ChN} mean poser in the range [f0-Bw/2, f0+Bw/2] 
  vector<string> varName; //ntuple variable naming
  varName.push_back("mode");
  varName.push_back("cycle");  
  varName.push_back("time");
  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
    stringstream tmp;
    tmp << iCh;
    varName.push_back("Ch"+tmp.str());
  }
  r_4 valTuple[varName.size()];
  NTuple powerEveloution(varName); 


  string fileName = "./" + sourceName_ + "/calibMonitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";

  POutPersist calibMonitorFile(fileName);
  
  //-----------------
  //Start real process
  //------------------
  for (int 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 = "./" + 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;
	valTuple[index++] = iMode;
	valTuple[index++] = icycle;
	valTuple[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();
	  valTuple[index++] = mean;
	}//end of channel loop
       
	//Fill Tuple
	powerEveloution.Fill(valTuple);
      }//end of files loop

    }//end of cycles loop
  }//end of mode loop

  //store Monitoring Tuple
  calibMonitorFile << PPFNameTag("monitor") <<  powerEveloution;

  //Compute Calibration Coefficient 


  return rc;
}//ProcessCalibration::processCmd

