//  
//  deplVandLumiExtract.C
//  macro to extract the depletion voltage measurements after production
//  and the luminosity information from the raw data
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <fstream>
#include <string.h>
#include <vector>
#include <map>
#include <locale>


//  ROOT include
#include "../../include/incBasicROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incHistROOT.h"
#include "../../include/incIOROOT.h"

#include "STDatabase.h"


#include "../../include/basicROOTTools.h"
#include "../../include/basicRooFitTools.h"

#include "../../include/lhcbstyle.C"








int deplVandLumiExtract(int, char**);

int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication* theApp = new TApplication("Analysis", &argc, argv,0,-1);
  argc = theApp->Argc();
  argv = theApp->Argv();

  bool runBatch = false;
  
  for (int i = 1; i<argc; i++) {
    if (strcmp(argv[i],"-b")==0) {
      runBatch = true;
      gROOT->SetBatch();
    }
  }
  
  style->SetPadLeftMargin(0.12);
  style->SetPadRightMargin(0.22);
  style->SetPadBottomMargin(0.15);
  style->SetPadTopMargin(0.05);
  
  int exit = deplVandLumiExtract(argc,argv);
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  return exit;
}


//  Executable method
int deplVandLumiExtract(int argc, char* argv[]){
  
  bool useRadius = false;
  
  TString det("TT");
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }
  }

  
  
  
  TString lay;
  if (strcmp(det.Data(),"TT")==0) {
    lay = "TTaU";
  }else{
    lay = "T3X2";
  }
  
  
  Info("deplVandLumiExtrac","Extract V_depl for %s / %s",lay.Data(),det.Data());
  
  
  
  const char* homeVar = std::getenv ("CCEHOME");
  const char* diskVar = std::getenv ("DISK");
  

  TString dn_lumi(Form("%s/data/lumi",homeVar));
  TString fn_lumi(Form("%s/lumi.root",
                       dn_lumi.Data()));
  
  TString dn_out(Form("%s/data/db",homeVar));
  TString fn_vdepl(Form("%s/vdepl.root",
                        dn_out.Data()));
  TString fn_flux(Form("%s/lumi.root",
                       dn_out.Data()));
  
  TString fn_db(Form("%s/%s_db.txt",
                      dn_out.Data(),
                      det.Data()));
  
  Printf("%s",fn_db.Data());
  
  std::ifstream f_db(fn_db.Data());
  
  std::vector<double> v_sector;
  std::vector<double> v_vdepl;
  
  std::vector<double> v_lumi7;
  std::vector<double> v_lumi8;
  std::vector<double> v_lumi13;
  std::vector<double> v_time;
  std::vector<double> v_fill;
  
  while (f_db) {
    std::string line;
    if (!getline(f_db,line)) {
      break;
    }
    std::istringstream s_line(line);
    std::vector<std::string> v_line;
    while (s_line) {
      std::string bit;
      if (!getline(s_line,bit,'\t')) {
        break;
      }
      v_line.push_back(bit);
    }
    if (v_line.size()>6 && strcmp(v_line[4].c_str(),lay.Data())==0) {
      Info("deplVandLumiExtract","Using sector %s",v_line[0].c_str());
      v_sector.push_back(atof(v_line[0].c_str()));
      v_vdepl.push_back(atof(v_line[6].c_str()));
    }else{
      Warning("deplVandLumiExtract","Database entry of %s has too few entries",v_line[0].c_str());
    }
  }
  
  for (int i=0; i<v_sector.size(); i++) {
    Info("deplVandLumiExtract","Sector %d has %6.2f V",
         (int)v_sector[i],v_vdepl[i]);
  }
  
  TFile* f_vdepl = new TFile(fn_vdepl.Data(),"UPDATE");
  
  TVectorD vr_sector  = getROOTVector<double>(v_sector);
  TVectorD vr_vdepl   = getROOTVector<double>(v_vdepl);
  
  f_vdepl->Delete(Form("v_sector_%s;*", lay.Data()));
  f_vdepl->Delete(Form("v_vdepl_%s;*",lay.Data()));
  
  vr_sector.Write(Form("v_sector_%s", lay.Data()));
  vr_vdepl.Write(Form("v_vdepl_%s",lay.Data()));

  
  f_vdepl->Close();
  
  Info("deplVandLumiExtract","Results written to %s",fn_vdepl.Data());
  
  
  TFile* f_lumi = TFile::Open(fn_lumi.Data());
  
  TVectorF vr_fillTime   = *(TVectorF*)f_lumi->Get("v_fillTime_tot");
  TVectorF vr_stopTime   = *(TVectorF*)f_lumi->Get("v_stopTime_tot");
  TVectorF vr_beamEner   = *(TVectorF*)f_lumi->Get("v_beamEner_tot");
  TVectorF vr_peakLumi   = *(TVectorF*)f_lumi->Get("v_peakLumi_tot");
  
  
  
  if (vr_fillTime.GetNoElements()!=vr_beamEner.GetNoElements() ||
      vr_fillTime.GetNoElements()!=vr_peakLumi.GetNoElements()) {
    Error("deplVandLumiExtract","Luminosity vectors do not have the same length");
    return EXIT_FAILURE;
  }
  
  double lumi7 = 0.0;
  double lumi8 = 0.0;
  double lumi13 = 0.0;
  
  
  
  std::map<int,double> timeMap = STDatabase::runTimeMap;
  
  for (int i=0; i<vr_fillTime.GetNoElements(); i++) {
    time_t raw_time = (long)vr_fillTime(i);
    if (vr_beamEner(i)==7.0) {
      lumi7 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i);
    }else if(vr_beamEner(i)==8.0){
      lumi8 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i);
    }else if(vr_beamEner(i)==13.0 || vr_beamEner(i)==14.0){
      lumi13 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i);
    }else{
      Warning("deplVandLumiExtract","Unknown center of mass energy for %s",ctime(&raw_time));
    }
    
    std::map<int,double>::iterator iter = timeMap.begin();
    for (; iter!=timeMap.end(); ++iter) {
      if (vr_fillTime(i)>iter->second) {
        v_fill.push_back((double)iter->first);
        v_lumi7.push_back(lumi7);
        v_lumi8.push_back(lumi8);
        v_lumi13.push_back(lumi13);
        Info("deplVandLumiExtract","At fill %d %s",iter->first,ctime(&raw_time));
        Printf("                   Luminosity @  7 TeV: %7.2f fb^{-1}",lumi7/1e9);
        Printf("                   Luminosity @  8 TeV: %7.2f fb^{-1}",lumi8/1e9);
        Printf("                   Luminosity @  13 TeV: %7.2f fb^{-1}",lumi13/1e9);
        v_time.push_back((double)iter->second);
        timeMap.erase(iter);
      }
    }
  }
  
  std::map<int,double>::iterator iter = timeMap.begin();
  for (; iter!=timeMap.end(); ++iter) {
    time_t fill_time = (long)iter->second;
    v_fill.push_back((double)iter->first);
    v_lumi7.push_back(lumi7);
    v_lumi8.push_back(lumi8);
    v_lumi13.push_back(lumi13);
    Info("deplVandLumiExtract","At fill %d: %s",iter->first,ctime(&fill_time));
    Printf("                   Luminosity @  7 TeV: %7.2f fb^{-1}",lumi7/1e9);
    Printf("                   Luminosity @  8 TeV: %7.2f fb^{-1}",lumi8/1e9);
    Printf("                   Luminosity @  13 TeV: %7.2f fb^{-1}",lumi13/1e9);
    v_time.push_back((double)iter->second);
    timeMap.erase(iter);
  }
  
  TFile* f_flux = new TFile(fn_flux.Data(),"UPDATE");
  
  TVectorD vr_time    = getROOTVector<double>(v_time);
  TVectorD vr_lumi7   = getROOTVector<double>(v_lumi7);
  TVectorD vr_lumi8   = getROOTVector<double>(v_lumi8);
  TVectorD vr_lumi13   = getROOTVector<double>(v_lumi13);
  TVectorD vr_fill    = getROOTVector<double>(v_fill);
  
  f_flux->Delete("v_time;*");
  f_flux->Delete("v_lumi7;*");
  f_flux->Delete("v_lumi8;*");
  f_flux->Delete("v_lumi13;*");
  f_flux->Delete("v_fill;*");
  
  vr_time.Write("v_time");
  vr_lumi7.Write("v_lumi7");
  vr_lumi8.Write("v_lumi8");
  vr_lumi13.Write("v_lumi13");
  vr_fill.Write("v_fill");
  
  f_flux->Close();
  
  Info("deplVandLumiExtract","Results written to %s",fn_flux.Data());
  
  
  return EXIT_SUCCESS;
}