/** To compute the sensitivity at a given FPS gs1:
 * - > Cd /opt/w3/MonitoringWeb/General/ 
 * - > vega -b
 * 
 * OR to copy the code to your own location dir: 
 * > cd dir
 * > cp  /opt/w3/MonitoringWeb/General/Sensitivity.C .
 * > cp  /opt/w3/MonitoringWeb/General/vegalogon.C .
 * > vega -b
 *
 * --------------------------------------------------------------------
 * Once VEGA is started, compile and start the script:
 * > .L Sensitivity.C++  (this will compile the Sensitivity macro)
 * > Sensitivity(gps1+300,300,....)   (one should change the 'path' argument)
 * 
 * --------------------------------------------------------------------
 * 
 */


#if !defined (__CINT__)

//includes
#include <iomanip>
#include <iostream>
#include <fstream>
#include <string>
#include <glob.h>
#include <unistd.h> // for the sleep function
#include <string>

//ROOT includes
#include <TROOT.h>
#include <TSystem.h>
#include <TApplication.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TLegend.h>
#include <TText.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TGraph.h>
#include <TImage.h>
#include <TImageDump.h>

//Virgo includes
#include <FrameL.h>
#include <Frv.h>

#endif

typedef enum {Virgo,VirgoP,VirgoP_withMS,AdV} VirgoConfigs;

void Sensitivity(//Time window definition
		 double gps = 0,        //end gps of the data to read (if zero, last stored gps)
		 double tAverage = 300, //duration of the data to read, in seconds
		 //FFT definition
		 double lengthFFT = 10,
		 //Quality criteria
		 int minStep = 12,
		 //Name of the server to be started
		 const char* serverPath = "/virgoLog/Cali",
		 const char* serverName = "TrendSensitivity"
		 );
// 		 bool saveImages = true,
// 		 //New clean data
// 		 const char* rawFileName   = 
// 		 const char* trendFileName = "/virgoData/ffl/trend.ffl",
// 		 //Output parameters for images
// 		 const char* path = "/opt/w3/MonitoringWeb/General",
// 		 //const char* path = "/olusers/virgorun",
// 		 //Base name for plots
// 		 const char* fileNameBase = "Sensitivity",
// 		 const char* historyFileNameBase = "General_Sensitivity",
// 		 //Feference sensitivity for plots
// 		 const char* referenceName = "July 12th, 2009",
// 		 const char* referenceFileName = "/virgoData/Cali/2009/2009_07/Sensitivity_931464610_2009Jul12_22h10m00.root"
// 		 );
// 		 //Plot range parameters
// 		 double fmin = 10.,
// 		 double fMax = 1e4,
// 		 double yMin = 1.e-23,
// 		 double yMax = 0.33e-18);

void End();
void print_log(//std::ofstream& logfile,
	       const char* line,
	       bool closefile=false);
int write_cfgfile(double tStart,
		  const char* serverPath,
		  const char* cfgfilename,	  
		  double tAverage,
		  double lengthFFT,
		  int minStep,
		  double gps);
const char* GetTFFileName(double tStart, double& gps_TF, double& gps_Sens);
std::string GetTimeChar(double tStart);

void CompareSensitivities(double gps, double tAverage, double lengthFFT,
			  const char* sensitivityROOTfile);
			  //			  bool saveImages,
// 			  double fmin,
// 			  double fMax,
// 			  double yMin,
// 			  double yMax,
// 			  const char* path = "/opt/w3/MonitoringWeb/Reconstruction",
// 			  const char* fileNameBase = "Sensitivities",
// 			  const char* historyFileNameBase = "Reconstruction_Sensitivities");

#if !defined (__CINT__)


// global variables
VirgoConfigs VirgoConf = VirgoP; //define the design sensitivity to be plotted

const char* CaliPath = "/virgoDev/Cali/v0r3p6"; 

const char* TiMPath  = "/virgo/VCS-7.0/VIRGOSW/TiM/v5r11p5";

//Input data files
const char* rawFileName   = "/virgoData/ffl/raw.ffl";
const char* trendFileName = "/virgoData/ffl/trend.ffl";
const char* hrecFileName  = "/virgoData/ffl/raw.ffl";

//Output paths and names
// Sensitivity 
const char* path = "/opt/w3/MonitoringWeb/General"; //base directory
const char* fileNameBase = "Sensitivity"; //Base name for plots
const char* historyFileNameBase = "General_Sensitivity"; //Base name for history plots
// Comparison with hrec sensitivity 
const char* pathHrec = "/opt/w3/MonitoringWeb/Reconstruction";//base directory
const char* fileNameBaseHrec = "Sensitivities"; //Banse name for plots
const char* historyFileNameBaseHrec = "Reconstruction_Sensitivities"; //Base name for history plots


//Plot range
double fMin = 10.;
double fMax = 10000.;
double yMin = 1.e-23;
double yMax = 0.33e-18;

//Reference sensitivity for plots
const char* referenceName = "July 12th, 2009";
const char* referenceFileName = "/virgoData/Cali/2009/2009_07/Sensitivity_931464610_2009Jul12_22h10m00.root";


bool saveImages = true;
bool savelog = true;


std::ofstream logfile; //log file for online checks


/** Generate a plot with the last sensitivity curve
 *  superpose onto a reference sensitivity curve 
 * 
 * \author Loic Rolland
 * 
 * \param gps end of the science mode data to use to compute sensitivity (if zero, last GPS stored in the data file)
 * \param tAverage duration in seconds of the data to be read to compute the averaged FFT
 * \param lengthFFT time window in seconds of a single FFT
 * 
 * \param minStep minimum step of the data to be analysed
 * 
 * \param serverName name of the server to be started
 *
 * \param saveImages do we save the gif images in the directory \param path ?
 * \param savelog do we save a log file called Sensitivity.log in directory \param path ?
 *  
 * \param rawFileName name of the file containing the raw data (channel h_20kHz)
 * \param trendFielName name of the file containing the trend data
 *
 * \param path path to save the output plots
 * \param fileNameBase base name for the output files for web use
 * \param historyFileNameBase base name for the output history files
 * 
 * \param referenceName name given to the reference sensitivity in the legend
 * \param referenceFileName name of the file with the reference sensitivity in h-strain, in vect format 
 *        (ABSOLUTE PATH!)
 * 
 * \param fMin, \param fMax plot x-axis frequency range [Hz]
 * \param yMin, \param yMax minimum and maximum y-axis sensitivity value [h/sqrt(Hz)]
 *
 */
void Sensitivity(//Time window definition
		 double gps,        //end gps of the data to read (if zero, last stored gps)
		 double tAverage, //duration of the data to read, in seconds
		 //FFT definition
		 double lengthFFT,
		 //Quality criteria
		 int minStep,
		 //Name of the server to be started
		 const char* serverPath,
		 const char* serverName)
// 		 bool saveImages,
// 		 bool savelog,
// 		 //New clean data
// 		 const char* rawFileName,
// 		 const char* trendFileName,
// 		 //Output parameters for images
// 		 const char* path,
// 		 //Base name for plots
// 		 const char* fileNameBase,
// 		 const char* historyFileNameBase,
// 		 //Feference sensitivity for plots
// 		 const char* referenceName,
// 		 const char* referenceFileName,
// 		 //Plot range paramters
// 		 double fMin,
// 		 double fMax,
// 		 double yMin,
// 		 double yMax)
{
  char fname[200];  
  char line[200];
  
  if(savelog) {
    char logfilename[200];
    sprintf(logfilename,"%s/Sensitivity.log",path);
    std::cout<<"Writing log file: "<<logfilename<<std::endl;
    logfile.open(logfilename);
  } 
  
  //window to search for data back in time (40 min)
  double tWindow = 40.*60; //s
  
  gROOT->SetStyle("Plain");
  gStyle->SetOptFit(111);
  gStyle->SetPalette(1);
  
  sprintf(line,"Sensitivity.C - Looking for good quality data up to GPS %9.0f\n",gps);
  print_log(line);
  print_log("Loading raw data and trend data");
  
  //Open the raw data file and get its time window
  FrFile *rawFile = FrFileINew(const_cast<char*>(rawFileName));
  double rawFile_tStart = FrFileITStart(rawFile);
  double rawFile_tEnd   = FrFileITEnd(rawFile);
  
  //Open the trend data file and get its time window
  FrFile *trendFile = FrFileINew(const_cast<char*>(trendFileName));
  double trendFile_tStart = FrFileITStart(trendFile);
  double trendFile_tEnd = FrFileITEnd(trendFile);    
  
  print_log("Searching for good quality data");
  sprintf(line,"Searching for Step %d data",minStep);
  print_log(line);
  //Define the end of the sensitivity clean data
  double tEnd = 0; //seconds
  if(gps == 0.) // default end time
    tEnd = rawFile_tEnd<trendFile_tEnd ? rawFile_tEnd : trendFile_tEnd ;
  else { //user defined end time
    tEnd = gps;
    if(tEnd > rawFile_tEnd || tEnd > trendFile_tEnd) {
      print_log("Sensitivity.C - ERROR: data for requested GPS time not available yet",true);
      End();return;
    }
  }
  
  //Define start time
  double tStart    = tEnd - tAverage; //later possible start time
  double tStartRef = tStart - tWindow; //sooner possible start time
  double tEndRef   = tEnd;
  
  //Check Reconstruction Quality flag
  FrVect *vAlpStatus = FrFileIGetV(trendFile,"Alp_Main_LOCK_STEP_STATUS", 
				   tStart, tAverage);
  
  double minstatus, maxstatus;
  FrVectMinMax(vAlpStatus, &minstatus, &maxstatus);
  FrVectFree(vAlpStatus);
  
  //Go back in time until data quality is correct,
  // until good quality data are found,
  // or the tAverage seconds window is no in the tWindow seconds window any more->no good quality data
  // search a window such that its end is 3 minutes before the unlock
  double tLarge = tAverage+180. ; //seconds
  while(minstatus < minStep) {// || meanOG < 7.e-10) {
    tEnd = tEnd - tAverage/10.;
    tStart = tEnd - tLarge;
    if(tEnd < tStartRef || tStart < rawFile_tStart || tStart < trendFile_tStart) {
      sprintf(line,"Sensitivity.C - ERROR: could not find enough good quality data in Step %d between GPS %9.0f and %9.0f\n",
	      minStep,tStart,tEndRef);
      print_log(line,true);
      End(); return;
    }
    
    FrVect *vAlpStatus = FrFileIGetV(trendFile,"Alp_Main_LOCK_STEP_STATUS", 
				     tStart, tLarge);
    FrVectMinMax(vAlpStatus, &minstatus, &maxstatus);
    FrVectFree(vAlpStatus);
  }
  sprintf(line,"Found good quality data at GPS %9.0f for %.0f s",tStart,tAverage);
  print_log(line);
  
  FrFileIEnd(rawFile);
  FrFileIEnd(trendFile);
  
  
  // //////////////////////////////////////////////////////
  //  COMPUTE THE SENSITIVITY IN THE CORRECT TIME WINDOW
  // //////////////////////////////////////////////////////
  
  //  const char* serverName = "OnlineSensitivity";   // name of the server to be sarted
  char command[200];  
  
  // calling the usual calibration executable (same as for online horizon)
  // -> write the cfg file with the correct GPS
  char cfgfilename[100];
  sprintf(cfgfilename,"%s.cfg",serverName);
  if(write_cfgfile(tStart,
		   serverPath,cfgfilename,
		   tAverage,lengthFFT,
		   minStep,
		   gps)!=0 ) return;

  {
    //Remove the output file if necessary
    sprintf(command,"rm -f %s/%s/HorizonTF_Sensitivity.root",serverPath,serverName);
    print_log(command);
    system(command);
    // -> start the process
    sprintf(command,"%s/${UNAME}/VirgoHorizonOnline.exe %s ",CaliPath,cfgfilename);
    print_log(command);
    system(command);
  }
  
  //wait until the output file with sensitivity is computed
  // (time out after 5 minutes)
  sprintf(fname,"%s/%s/HorizonTF_Sensitivity.root",serverPath,serverName);
  int timeout = 10*60; //seconds
  sprintf(line,"Looking for the output sensitivity curve %s (with timeout of %f min)",
	  fname,timeout/60.);
  print_log(line);
  int time = 0;
  int dtime = 15;     //seconds
  TFile* file = 0;
  while(time<timeout && file==NULL) {
    sleep(dtime);
    time += dtime;
    sprintf(line,"  ... Waiting for the ROOT file since %d s",time);
    print_log(line);
    delete file;
    file = new TFile(fname);
    gROOT->cd();
    if(!file || !file->IsOpen()) {
      delete file;
      file = NULL;
    }
  }
  if(!file || !file->IsOpen()) {
    sprintf(line,"Could not find ROOT file with sensitivity curve: %s (after timeout %f min)",fname,timeout/60.);
    print_log(line,true);
    End();return;
  }
  std::cout<<"Found ROOT sensitivity file: "<<file<<std::endl;
  print_log("Found ROOT sensitivity file");
  TH1D* hsensitivity = (TH1D*)file->Get("FFT <amplitude>(V1:Pr_B1_ACp)"); //old name
  if(!hsensitivity)
    hsensitivity = (TH1D*)file->Get("hSensitivityH");
  file->Clear();
  file->Close();
  delete file;
  
  if(!hsensitivity) {
    sprintf(line,"Could not find the sensitivity curve in ROOT file !");
    print_log( line,true);
    End();return;
  }
  

  // ////////////////////////////////////////////////////////////////////
  // PLOT THE CURRENT SENSITIVITY CURVE AND THE REFERENCE ONE
  // ////////////////////////////////////////////////////////////////////
  TCanvas *can = new TCanvas("Sensitivity","Sensitivity",10,10,900,700);
  // Plot sensitivity 
  can->SetBorderMode(0);
  can->SetFillColor(0);
  can->SetLogx(1);
  can->SetLogy(1);
  can->SetGridx(1);
  can->SetGridy(1);
  gStyle->SetOptStat(0);
  char drawOpt[20];
  sprintf(drawOpt,"SAME");
  
  //Plot range (empty) histogram
  TH2F* hrange = new TH2F("hrange","hrange",100,fMin,fMax,100,yMin,yMax);
  hrange->SetBit(TH1::kNoStats);
  hrange->SetBit(TH1::kNoTitle);
  hrange->GetXaxis()->SetTitle("Frequency [Hz]");
  hrange->GetYaxis()->SetTitle("h / #sqrt{Hz}");
  hrange->GetXaxis()->SetTitleOffset(1.1);
  hrange->GetYaxis()->SetTitleOffset(1.3);
  hrange->Draw();
  
  //Plot design sensitivity curve
  const char* designFileNameDir = "/opt/w3/DataAnalysis/Calibration/Sensitivity/Design";
  char designFileName[200];
  switch(VirgoConf) {
  case Virgo:
    sprintf(designFileName,"%s/VIRGO_DesignSensitivityH_nolines.txt",designFileNameDir);
    break;
  case VirgoP:
    sprintf(designFileName,"%s/VIRGOPLUS_DesignSensitivityH.txt",designFileNameDir);
    break;
  case VirgoP_withMS:
    sprintf(designFileName,"%s/VIRGOPLUS_MS_DesignSensitivityH.txt",designFileNameDir);
    break;
  case AdV: //not yet ready
    sprintf(designFileName,"%s/ADVVIRGO_DesignSensitivityH.txt",designFileNameDir);
    break;
  }
  sprintf(line,"Looking for the design sensitivity curve from %s",designFileName);
  if(savelog) print_log(line);
  TGraph* gDesign = new TGraph(designFileName);
  if(gDesign) {
    std::cout<<"Drawing design "<<gDesign<<std::endl;
    gDesign->SetLineColor(15);
    gDesign->SetMarkerColor(15);
    gDesign->Draw("LSAME");
  }
  else {
    sprintf(line,"Could not build the TGraph design sensitivity curve from file %s",designFileName);
    if(savelog) print_log(line);
  }
  
  //Plot reference sensitivity
  sprintf(line,"Adding the reference sensitivity curved from %s",referenceFileName);
  print_log(line);
  TH1D* refSens = NULL;
  TFile *fileRef = new TFile(referenceFileName);
  gROOT->cd();
  if(!fileRef) {
    sprintf(line,"Could not open the file %s",referenceFileName);
    print_log(line);    
  }
  else {
    TH1D *refSensTmp = (TH1D *)fileRef->Get("hSensitivityH");
    if(!refSensTmp) refSensTmp = (TH1D*)fileRef->Get("Sensitivity");
    fileRef->Clear();
    fileRef->Close();
    delete fileRef;
    if(!refSensTmp) {
      sprintf(line,"Could not find TH1D hSensitivityH in the file %s",referenceFileName);
      print_log(line);    
    }
    
    //Copy the VSPlot into a TH1D... else seg fault when saving the plot into an image !
    refSens = new TH1D(*refSensTmp);
    if(refSens) {
      std::cout<<"Drawing refsens "<<refSens<<std::endl;
      refSens->SetBit(TH1::kNoStats);
      refSens->SetBit(TH1::kNoTitle);
      refSens->GetXaxis()->SetRangeUser(fMin,fMax);
      refSens->GetXaxis()->SetTitle("Frequency [Hz]");
      refSens->GetYaxis()->SetTitle("h / #sqrt{Hz}");
      refSens->SetMinimum(yMin);
      refSens->SetMaximum(yMax);
      refSens->Draw(drawOpt);
      sprintf(drawOpt,"SAME");
    }
  }
  
  //Plot the current sensitivity
  std::cout<<"Drawing the hsensitivity "<<hsensitivity<<std::endl;
  
  if(hsensitivity) {
    hsensitivity->GetXaxis()->SetRangeUser(fMin,fMax);
    hsensitivity->GetXaxis()->SetTitle("Frequency [Hz]");
    hsensitivity->GetYaxis()->SetTitle("h / #sqrt{Hz}");
    hsensitivity->SetMinimum(yMin);
    hsensitivity->SetMaximum(yMax);  
    hsensitivity->SetBit(TH1::kNoStats);
    hsensitivity->SetLineColor(2);
    hsensitivity->SetLineWidth(2);
    hsensitivity->Draw(drawOpt);
    sprintf(drawOpt,"SAME");
  }
  
  //Get the time into a char
  std::string str = GetTimeChar(tStart);
  char sensTime[150];
  sprintf(sensTime,"%s",str.c_str());
  str.clear();
  
  // Prepare for archiving  
  char refSensTime[200];
  sprintf(refSensTime,"Reference Sensitivity: %s",referenceName);
  
  TLegend* leg1 = new TLegend(.28,.78,.88,.90); 
  leg1->Draw();
  leg1->SetFillStyle(1001);
  leg1->SetFillColor(10);
  if(hsensitivity) leg1->AddEntry(hsensitivity,sensTime,"L");
  if(refSens)      leg1->AddEntry(refSens,refSensTime,"L");  
  char designLeg[50];
  switch(VirgoConf) {
  case Virgo:
    sprintf(designLeg,"Virgo design");
    break;
  case VirgoP:
    sprintf(designLeg,"Virgo+ design");
    break;
  case VirgoP_withMS:
    sprintf(designLeg,"Virgo+ with MS design");
    break;
  case AdV: //not yet ready
    sprintf(designLeg,"Adv Virgo design");
    break;
  }
  if(gDesign)  leg1->AddEntry(gDesign,designLeg,"L");
  
  //   TText* text1 = new TText(100,1e-19,"NOT YET CALIBRATED !");
  //   text1->Draw();
  
  can->Modified();
  can->Update();
  
  if(saveImages) { //save plots into files
    sprintf(line,"Saving images in %s",path);
    print_log(line);
    
    //Create history directory
    sprintf(command,"mkdir %s/history",path);
    system(command);
    
    //Get the date for file names
    time_t tt = time_t(tStart + FRGPSOFF);
    char date[40];
    strftime(date,24,"%d-%b-%Y-%Hh%M",localtime(&tt));
    
    //Save the jpg image
    print_log("Saving the jpg image");
    char jpgFileName[150];
    sprintf(jpgFileName,"%s/plots/%s.jpg",path, fileNameBase);
    TImageDump *imgdump = new TImageDump(jpgFileName); 
    print_log("Painting the canvas");
    can->Paint(); 
    print_log("Done");
    TImage *img = imgdump->GetImage(); // get internal image 
    img->SetImageQuality(TAttImage::kImgGood); 
    img->SetImageCompression(60); 
    imgdump->Close(); // that will delete img
    print_log("Saved done");
    
    //Copy the jpg image in archive directory
    char historyJpgFileName[150];
    sprintf(historyJpgFileName,"%s/history/%s-%9.0f-%s.jpg",path, historyFileNameBase, tStart, date);
    gSystem->Exec(Form("cp %s %s", jpgFileName, historyJpgFileName));
    
    //Save the small jpg image
    char smallJpgFileName[150];
    char execConvertSmall[500];
    sprintf(smallJpgFileName,"%s/plots/%s_small.jpg",path, fileNameBase);
    sprintf(execConvertSmall,"convert -geometry 386x300 -density 72x72 %s %s",jpgFileName, smallJpgFileName);
    gSystem->Exec(execConvertSmall);
    
    //Save in ascii file
    sprintf(historyJpgFileName,"%s/history/%s-%9.0f-%s.txt",path, historyFileNameBase, tStart, date);
    sprintf(line,"Writing ascii file %s",historyJpgFileName);
    print_log(line);
    std::ofstream txtfile(historyJpgFileName);
    if(txtfile) {
      for(int bin = 1 ; bin<hsensitivity->GetNbinsX(); bin++) {
	double freq = hsensitivity->GetBinCenter(bin);
	if(freq>=5 && freq<=10000)
	  txtfile<< freq <<"  "<<hsensitivity->GetBinContent(bin)<<std::endl;
      }
      txtfile.flush();
      txtfile.close();
    }
  }
  
  print_log("Sensitivity.C - End, now starting CompareSensitivities");
  
  sprintf(fname,"%s/%s/HorizonTF_Sensitivity.root",serverPath,serverName);  
  CompareSensitivities(tStart,tAverage,lengthFFT,
		       fname);
// 		       saveImages,rawFileName,fMin,fMax,yMin,yMax);
  
  print_log("END",true);
  End();
}
  
void End()
{
  if(gROOT->IsBatch())
    gROOT->GetApplication()->Terminate(0);  
}

void print_log(const char* line,
	       bool closefile)
{
  std::cout<<line<<std::endl;
  if(!logfile || logfile.bad()) return;
  logfile<<line<<std::endl;
  if(closefile) {
    logfile.flush();
    logfile.close();
  }
}


/** Write the configuration file (.cfg) used as argument for the executable
 * computing the Virgo sensitivity and horizon
 * 
 * \param tStart start time of the data to compute the sensitivity
 * \param serverPath path of the server to start the log file (and output files)
 * \param cfgfilename name of the configuration file to write
 * \param tAverage duration of the data to average, in seconds
 * \param lengthFFT duration of a single FFT, in seconds
 * \param minStep minimum locking step to use the data
 * \param gps determine what Virgo TF will be used:
 *        \li 0:  use the TF stored in /virgoData/Cali/lastTF.vect
 *        \li >0: search for the last TF computed before the tStart GPS
 */
int write_cfgfile(double tStart,
		  const char* serverPath,
		  const char* cfgfilename,	  
		  double tAverage,
		  double lengthFFT,
		  int minStep,
		  double gps)
{
  std::ofstream ofile(cfgfilename);
  if(!ofile) return -1;
  
  ofile<<"CFG_PRIO 4               # Main priority; 0 means no change (nice(0))"<<std::endl;
  ofile<<"CFG_NOFILESAVE"<<std::endl;
  ofile<<"CFG_NODBSAVE              # No commit into Db"<<std::endl;
  ofile<<"CFG_PWD "<<serverPath<<"  # Current logfile path <path>/<cmName>"<<std::endl;
  ofile<<"CFG_RFLAG eloff           # Report flags: el{on,off}, stdout{on,off}, log{on,off}"<<std::endl;
  
  ofile<<"FDIN_FILE   /virgoData/ffl/raw.ffl  "<<std::setprecision(9)<<tStart<<" "<<tAverage<<std::endl;
  //ofile<<"FDIN_FRAME_DURATION 1"<<std::endl;
  
  ofile<<"HORIZON_DFCHANNEL       V1:Pr_B1_ACp"<<std::endl;
  if(gps<=0) 
    ofile<<"HORIZON_TFFILENAME       /virgoData/Cali/lastTF.vect"<<std::endl;
  else {
    double gps_TF=-1, gps_Sens=-1;
    const char* tffilename = GetTFFileName(tStart, gps_TF, gps_Sens);
    if(tffilename==NULL || gps_TF==-1 || gps_Sens==-1) {
      std::cerr<<"ERROR: could not find the correct TF file for tStart: "<<std::setprecision(9)<<tStart<<std::endl;
      std::cout<<tffilename<<" "<<gps_TF<<" "<<gps_Sens<<std::endl;
      return -1;
    }
    ofile<<"HORIZON_TFFILENAME      "<<tffilename<<std::endl;
    ofile<<"HORIZON_TF_GPS          "<<std::setprecision(9)<<gps_TF  <<std::endl;
    ofile<<"HORIZON_Sensitivity_GPS "<<std::setprecision(9)<<gps_Sens<<std::endl;    
  }
  ofile<<"HORIZON_FreqRefLine      379.0"<<std::endl; 
  
  ofile<<"HORIZON_DarmChannelSMS    V1:Gc"<<std::endl;
  ofile<<"HORIZON_DarmChannelParam  Sensing_Sine_Darm_Amp"<<std::endl;
  
  ofile<<"HORIZON_FFTLENGTH     "<< lengthFFT <<std::endl;
  ofile<<"HORIZON_AVERAGINGTIME "<< tAverage  <<std::endl;
  ofile<<"HORIZON_VERBOSELEVEL  1"<<std::endl;
  
  ofile<<"HORIZON_MINSTEP       "<<minStep<<std::endl;
  
  ofile<<"HORIZON_SaveROOTFile true"<<std::endl;
  
  return 0;
}

/** Get the GPS times of the last TF and Virgo Sensitivity measurements
 *
 */
const char* GetTFFileName(double tStart, double& gps_TF, double& gps_Sens)
{ 
  //Get the last TF and sensitivity GPS in the log file
  std::ifstream file("/virgoData/Cali/ActuatorGainCalibrationLogs/TFSensitivityData.txt");
  if(file && file.good()) {
    double gpsTFStart, gpsTFEnd, gpsSensStart, gpsSensEnd;
    double lastgpsTFStart=-1, lastgpsSensStart=-1;
    char tmpchar[20];
    while(1) {
      file>>tmpchar>>gpsTFStart>>gpsTFEnd>>gpsSensStart>>gpsSensEnd;
      if(!file.good()) break;
      if(gpsTFStart>tStart && gpsSensStart>tStart) break; //Found a TF/Sens made AFTER the interesting GPS
      //Store the GPS of the last TF/Sens if both were computed
      if(gpsTFStart>0 && gpsSensStart>0) {
	lastgpsTFStart = gpsTFStart;
	lastgpsSensStart = gpsSensStart;
      }
    }
    file.close();
    if(lastgpsTFStart>0 && lastgpsSensStart>0) {
      gps_TF   = lastgpsTFStart;
      gps_Sens = lastgpsSensStart;
      //Now search for the TF vect file
      glob_t globbuf;
      char path[200];
      sprintf(path,"/virgoData/Cali/200\?/200\?_\?\?/TF_%9.0f_*.vect",gps_TF);
      std::cout << "\033[34;1m Searching for " << path <<"\033[0m" << std::endl;      
      if(glob(path, GLOB_ERR, 0, &globbuf) == 0 && globbuf.gl_pathc != 0) {
	if(globbuf.gl_pathc>1) 
	  std::cout<<"WARNING: Found "<<globbuf.gl_pathc<<" files for the Virgo TF at"<<std::setprecision(9)<<gps_TF<<std::endl;
	std::cout<<"\033[35;1m -> Found TF file "<<globbuf.gl_pathv[0]<<"\033[0m"<<std::endl;
	return globbuf.gl_pathv[0];
      }
      else {
	std::cerr<<"ERROR: could not find Virgo TF file !!!"<<std::endl;
	return NULL;
      }
    }
  }
  return NULL;  
}



std::string GetTimeChar(double tStart)
{
  //Get the date with 'nice' format for the plot legend
  char gpsCommand[1000];
  sprintf(gpsCommand,"%s/$UNAME/TiMGpsTime.exe -gps %9.0f > /virgoData/Md/string.txt \n",
	  TiMPath,tStart);
  char line[200];
  sprintf(line,"Sensitivity.C - Executing %s\n",gpsCommand);
  print_log(line);
  
  gSystem->Exec(gpsCommand);
  
  FILE *f = fopen("/virgoData/Md/string.txt","r");
  char sensTime[200];
  char toto[64];
  fscanf(f,"%s \n",sensTime);
  fscanf(f,"%s",toto);
  strcat(sensTime," ");
  strcat(sensTime,toto);
  fscanf(f,"%s",toto);
  fscanf(f,"%s",toto);
  strcat(sensTime," ");
  strcat(sensTime,toto);
  fscanf(f,"%s",toto);
  strcat(sensTime," ");
  strcat(sensTime,toto);
  fscanf(f,"%s",toto);
  strcat(sensTime," ");
  strcat(sensTime,toto);
  fclose(f);
  gSystem->Exec("rm /virgoData/Md/string.txt");
  
  sprintf(line,"GetTimeChar: %s",sensTime);
  print_log(line);
  return std::string(sensTime);
}

void CompareSensitivities(double gps, double tAverage, double lengthFFT,
			  const char* sensitivityROOTfile)
// 			  bool saveImages,
// 			  const char* rawFileName,
// 			  double fMin,
// 			  double fMax,
// 			  double yMin,
// 			  double yMax,
// 			  const char* path,
// 			  const char* fileNameBase,
// 			  const char* historyFileNameBase)
{
  if(gps<=0 || tAverage<0 || lengthFFT<=0 || sensitivityROOTfile==NULL)
    return;
  char* signalName = "h_20000Hz";
  
  char fname[200];
  char line[200];
  
  sprintf(line,"Looking for the output sensitivity curve");
  print_log(line);
  int time = 0;
  int timeout = 3*60; //seconds
  int dtime = 10;     //seconds
  TFile* file = 0;
  while(time<timeout && file==NULL) {
    sprintf(line,"Searching for the ROOT file (%d)",time);
    print_log(line);
    delete file;
    file = new TFile(sensitivityROOTfile);
    gROOT->cd();
    if(!file || !file->IsOpen()) {
      delete file;
      file = NULL;
      sleep(dtime);
      time += dtime;
    }
  }
  std::cout<<"Found ROOT sensitivity file: "<<file<<std::endl;
  if(!file || !file->IsOpen()) {
    sprintf(line,"Could not find ROOT file with sensitivity curve: %s",fname);
    print_log(line);
  }
  TH1D* hsensitivity = (TH1D*)file->Get("FFT <amplitude>(V1:Pr_B1_ACp)"); //old name
  if(!hsensitivity)
    hsensitivity = (TH1D*)file->Get("hSensitivityH");
  file->Clear();
  file->Close();
  delete file;  
  if(!hsensitivity) {
    sprintf(line,"Could not find the sensitivity curve in ROOT file !");
    print_log(line);
    return;
  }
  
  sprintf(line,"Computing the sensitivity from h(t) at GPS %9.0f for %.0f s",gps,tAverage);
  print_log(line);
  
  //Compute the sensitivity from h(t)
  FrvRFFT* fft = NULL;
  //Loop on the data to compute the FFTs
  double start = gps;
  int cpt = 0;
  FrFile* ifile = FrFileINew(const_cast<char*>(hrecFileName));
  if(!ifile) {
    sprintf(line,"Could not find the file %s",hrecFileName);
    print_log(line);
    return;
  }
  while(start+lengthFFT <= gps+tAverage) {
    //Get the data vectors     
    FrVect* vSignal = FrFileIGetV(ifile,signalName,start,lengthFFT);
    if(!vSignal) {
      sprintf(line,"Could not get the FrVect %s in file %s at GPS %9.0f for %03.1f s",
	      signalName,rawFileName,start,lengthFFT);
      print_log(line);
      return;
    }
    //Compute FFT
    fft = FrvRFFTFor(fft,vSignal);      
    //Free memory
    FrVectFree(vSignal);
    
    //prepare next step
    start += lengthFFT;
    cpt++;
    
    if(lengthFFT>=10 || (lengthFFT<10 && !(cpt%10))) std::cout<<"."<<std::flush;
  }
  if(!fft) return;
  
  //Make the histogram from the FFT
  FrVect* vect = fft->amplitudeA;
  if(!vect) return;
  //FrVect characteristics
  unsigned nBins = vect->nData;  
  double xStart = vect->startX[0];
  double binwidth = vect->dx[0];  
  //Create histos
  //Spectrum of h(t)
  TH1D* histo = new TH1D(vect->name,vect->name,nBins,xStart,xStart+nBins*binwidth); 
  histo->SetBit(TH1::kNoStats);
  histo->SetXTitle(vect->unitX[0]);
  histo->SetYTitle(vect->unitY);  
  //Ratio FFT[h(t)] / h(f)
  TH1D* hratio = new TH1D("FFThoft_over_hoff","FFThoft_over_hoff",
			  nBins,xStart,xStart+nBins*binwidth); 
  hratio->SetBit(TH1::kNoStats);
  hratio->SetXTitle(vect->unitX[0]);
  hratio->SetYTitle("Ratio");
  //Fill histogram
  for(unsigned i = 0 ; i<nBins ; i++) {
    double value = 0;
    if(     vect->type == FR_VECT_8R) value = vect->dataD[i];  // double
    else if(vect->type == FR_VECT_4R) value = vect->dataF[i];  // float
    else {
      sprintf(line,"VectToTH1D: data type not supported (%d,%s)",vect->type,vect->name);
      print_log(line);
      delete histo;
      return;
    }
    histo->SetBinContent(i+1, value); //since TH1D bins are from 1 to nBins included
    
    //Compute ratio
    double freq = xStart + (i+0.5)*binwidth; //Hz
    double hoff = hsensitivity->GetBinContent(hsensitivity->FindBin(freq));
    if(hoff>0)
      hratio->SetBinContent(i+1, value / hoff);
  }
  
  // //////////////////////////////////
  //Make the plot of the two sensitivities
  TCanvas* can = new TCanvas("Sensitivities","Sensitivities",
			     10,10,900,700);
  can->cd();
  // Plot sensitivity 
  can->SetBorderMode(0);
  can->SetFillColor(0);
  can->SetLogx(1);
  can->SetLogy(1);
  can->SetGridx(1);
  can->SetGridy(1);
  gStyle->SetOptStat(0);
  char drawOpt[20];
  sprintf(drawOpt,"L");
  
  //Get the time in a char
  std::string str = GetTimeChar(gps);
  char sensTime[150];
  sprintf(sensTime,"%s",str.c_str());
  str.clear();
  
  //plot the frequency-domain sensitivity curve
  hsensitivity->SetBit(TH1::kNoStats);
  hsensitivity->SetBit(TH1::kNoTitle);
  //  hsensitivity->SetTitle(sensTime);
  hsensitivity->GetXaxis()->SetRangeUser(fMin,fMax);
  hsensitivity->GetXaxis()->SetTitle("Frequency [Hz]");
  hsensitivity->GetYaxis()->SetTitle("h / #sqrt{Hz}");
  hsensitivity->SetMinimum(yMin);
  hsensitivity->SetMaximum(yMax);
  hsensitivity->SetLineColor(2); //red, same as in General trend data plot
  hsensitivity->Draw(drawOpt);
  sprintf(drawOpt,"LSAME");
  
  //plot the time-domain sensitivity curve
  histo->SetLineColor(4);  //blue, for better visibility than purple as usual
  histo->Draw(drawOpt);    
  
  TLegend* leg1 = new TLegend(.28,.78,.88,.90); 
  std::cout<<"Setting legend header to "<<sensTime<<std::endl;
  sprintf(line,"Setting legend header to %s",sensTime);
  print_log(line);
  leg1->SetHeader(sensTime);  
  leg1->SetFillStyle(1001);
  leg1->SetFillColor(10);
  if(hsensitivity) leg1->AddEntry(hsensitivity,"Frequency domain","L");
  if(histo)        leg1->AddEntry(histo,"h(t) spectrum","L");  
  leg1->Draw();
  
  can->Modified();
  can->Update();
  
    
  // ///////////////
  //Make the plot of the ratio
  TCanvas* can2 = new TCanvas("Ratio_HofT_vs_HofF","Ratio_HofT_vs_HofF",
			      10,10,900,700);
  can2->cd();
  // Plot sensitivity 
  can2->SetBorderMode(0);
  can2->SetFillColor(0);
  can2->SetLogx(1);
  can2->SetGridx(1);
  can2->SetGridy(1);
  gStyle->SetOptStat(0);
  
  //plot the frequency-domain sensitivity curve
  hratio->SetBit(TH1::kNoStats);
  hratio->SetBit(TH1::kNoTitle);
  hratio->GetXaxis()->SetRangeUser(fMin,fMax);
  hratio->GetXaxis()->SetTitle("Frequency [Hz]");
  hratio->GetYaxis()->SetTitle("FFT[h(t)] / h(f)");
  hratio->SetMinimum(0.2);
  hratio->SetMaximum(2.0);
  hratio->SetLineColor(2); //red, same as in General trend data plot
  hratio->Draw("L");
  
  TLegend* leg2 = new TLegend(.28,.78,.88,.90); 
  std::cout<<"Setting legend header to "<<sensTime<<std::endl;
  sprintf(line,"Setting legend header to %s",sensTime);
  print_log(line);
  leg2->SetHeader(sensTime);  
  leg2->SetFillStyle(1001);
  leg2->SetFillColor(10);
  //if(hratio) leg2->AddEntry(hratio,"#frac{FFT[h(t)]}{h(f)}","L");
  leg2->Draw();
  
  can2->Modified();
  can2->Update();
  

  
  // //////////////////////////////////////  
  //Save the images
  if(saveImages) { //save plots into files
    sprintf(line,"Saving images in %s",pathHrec);
    print_log(line);
    
    //Create history directory
    char command[200];
    sprintf(command,"mkdir %s/history",pathHrec);
    system(command);
    
    //Get the date for file names
    time_t tt = time_t(gps + FRGPSOFF);
    char date[40];
    strftime(date,24,"%d-%b-%Y-%Hh%M",localtime(&tt));
    
    //First image
    //Save the jpg image
    std::cout<<"Saving the jpg image"<<std::endl;
    char jpgFileName[150];
    sprintf(jpgFileName,"%s/plots/%s.jpg",pathHrec, fileNameBaseHrec);
    TImageDump *imgdump = new TImageDump(jpgFileName); 
    print_log("Painting the canvas");
    can->Paint(); 
    print_log("Done");
    TImage *img = imgdump->GetImage(); // get internal image 
    img->SetImageQuality(TAttImage::kImgGood); 
    img->SetImageCompression(60); 
    imgdump->Close(); // that will delete img
    print_log("Saved done");
    
    //Copy the jpg image in archive directory
    char historyJpgFileName[150];
    sprintf(historyJpgFileName,"%s/history/%s-%9.0f-%s.jpg",pathHrec, historyFileNameBaseHrec, gps, date);
    gSystem->Exec(Form("cp %s %s", jpgFileName, historyJpgFileName));
    
    //Save the small jpg image
    char smallJpgFileName[150];
    char execConvertSmall[500];
    sprintf(smallJpgFileName,"%s/plots/%s_small.jpg",pathHrec, fileNameBase);
    sprintf(execConvertSmall,"convert -geometry 386x300 -density 72x72 %s %s",jpgFileName, smallJpgFileName);
    gSystem->Exec(execConvertSmall);       
    
    //Second image
    //Save the jpg image
    std::cout<<"Saving the jpg image (2)"<<std::endl;
    sprintf(jpgFileName,"%s/plots/%s_ratio.jpg",pathHrec, fileNameBase);
    TImageDump *imgdump2 = new TImageDump(jpgFileName); 
    print_log("Painting the canvas (2)");
    can2->Paint(); 
    print_log("Done");
    TImage *img2 = imgdump2->GetImage(); // get internal image 
    img2->SetImageQuality(TAttImage::kImgGood); 
    img2->SetImageCompression(60); 
    imgdump2->Close(); // that will delete img
    print_log("Saved done");
    
    //Copy the jpg image in archive directory
    sprintf(historyJpgFileName,"%s/history/%s_ratio-%9.0f-%s.jpg",pathHrec, historyFileNameBaseHrec, gps, date);
    gSystem->Exec(Form("cp %s %s", jpgFileName, historyJpgFileName));
    
    //Save the small jpg image
    sprintf(smallJpgFileName,"%s/plots/%s_ratio_small.jpg",pathHrec, fileNameBaseHrec);
    sprintf(execConvertSmall,"convert -geometry 386x300 -density 72x72 %s %s",jpgFileName, smallJpgFileName);
    gSystem->Exec(execConvertSmall); 
  }  
}

#endif