//  
//  plotSector.C
//  Plot measured depletion voltage and Hamburg model for a certain read-out sector
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Updated by Elena Graverini on 10.06.15
//  University of Zurich, elena.graverini@cern.ch
//
//  compile with:
//  g++ -o plotSector plotSector.C `root-config --glibs --cflags`
//
//  Copyright only for commercial use

//  General include
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <stdlib.h>
#include <limits>
#include <string>
#include <time.h>


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


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


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


#include "deplTool.h"
#include "hamburgTool.h"


// -d detector (TT/IT)
// -s read-out sector number
// -r consider only hit in a distance closer than r mm from the beam axis
// -v number of summed up strips
// -ns Do not save the value in the root file


int plotSector(int argc, char* argv[]);


int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication* theApp = new TApplication("Analysis", &argc, argv);
  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();
    }
  }
  
  
  int exit = plotSector(argc,argv);
  
  
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  
  return exit;
}


//  Executable method
int plotSector(int argc, char* argv[]){
  
  TVectorD v_res(5);
  
  gStyle->SetPadLeftMargin(0.20);
  gStyle->SetTitleOffset(1.5,"Y");
  gStyle->SetTitleOffset(1.0,"X");

  //const char* diskVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";//std::getenv ("DISK");
  const char* diskVar = std::getenv ("DISK");
  //const char* homeVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";//std::getenv ("HOME");
  const char* homeVar = std::getenv ("CCEHOME");
  

  TString det("TT");
  TString lay("TTaU");
  int val        = 3;
  int sector     = 2634;
  double radius  = -1000;
  bool notSave   = false;
  
  
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) {
      det = argv[++i];
    }else if (strcmp(argv[i],"-v")==0 && i+1<argc) {
      val = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
      sector = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
      radius = std::atof(argv[++i]);
    }else if (strcmp(argv[i],"-ns")==0) {
      notSave = true;
    }
  }
  
  
  if (det.EqualTo("IT")) {
    lay = "T3X2";
  }else{
    lay = "TTaU";
  }
  
  // Directory and file name for output graph
  TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d",
                        diskVar,det.Data(),lay.Data(),sector));
  TString fn_graph(Form("%s/time_evolv",
                        dn_graph.Data(),val));
  
  if (radius>0.0) {
    fn_graph += Form("_r%d",(int)radius);
  }
  
  fn_graph += ".pdf";
  
  
  // Directory and file name for measured Vdepl, initial depletion voltage, luminosity and Hamburg model output
  TString dn_data(Form("%s/data/ST/Aging",diskVar));
  TString fn_data(Form("%s/CCEScan.root",dn_data.Data()));
  
  TString dn_deplV(Form("%s/data/db",homeVar));
  TString fn_deplV(Form("%s/vdepl.root",dn_deplV.Data()));
  
  TString dn_fill(Form("%s/data/db",homeVar));
  TString fn_fill(Form("%s/lumi.root",dn_fill.Data()));
  
  TString dn_hamburg(Form("%s/data/ST/Aging",diskVar));
  TString fn_hamburg(Form("%s/Hamburg.root",dn_hamburg.Data()));

  
  TFile* f_data    = TFile::Open(fn_data.Data());
  TFile* f_deplV   = TFile::Open(fn_deplV.Data());
  TFile* f_fill    = TFile::Open(fn_fill.Data());
  TFile* f_hamburg = TFile::Open(fn_hamburg.Data());
  

  TVectorD* vp_fill  = (TVectorD*)f_fill->Get("v_fill");
  TVectorD* vp_time  = (TVectorD*)f_fill->Get("v_time");

  
  if (!vp_fill || !vp_time) {
    Error("plotSector","Luminosity vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_fill    = *vp_fill;
  TVectorD v_time    = *vp_time;
  
  TVectorD* vp_sector  = (TVectorD*)f_deplV->Get(Form("v_sector_%s",lay.Data()));
  TVectorD* vp_vdepl   = (TVectorD*)f_deplV->Get(Form("v_vdepl_%s",lay.Data()));
  
  if (!vp_sector || !vp_vdepl) {
    Error("plotSector","Depletion Voltage vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_sector  = *vp_sector;
  TVectorD v_vdepl   = *vp_vdepl;
  
  TString vn_output = TString::Format("v_Hamburg_%s_%d",
                                      lay.Data(),sector);
  
  if (radius>0.0) {
    vn_output += Form("_r%d",(int)radius);
  }
  
  TVectorD* vp_deltaV  = (TVectorD*)f_hamburg->Get(Form("%s_deltaV",vn_output.Data()));
  TVectorD* vp_time_H  = (TVectorD*)f_hamburg->Get(Form("%s_time",  vn_output.Data()));
  TVectorD* vp_hig     = (TVectorD*)f_hamburg->Get(Form("%s_hig",   vn_output.Data()));
  TVectorD* vp_low     = (TVectorD*)f_hamburg->Get(Form("%s_low",   vn_output.Data()));
  
  if (!vp_deltaV || !vp_time_H || !vp_hig || !vp_low) {
    Error("plotSector","Hamburg model vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_deltaV   = *vp_deltaV;
  TVectorD v_time_H   = *vp_time_H;
  TVectorD v_hig      = *vp_hig;
  TVectorD v_low      = *vp_low;
  
  double vdeplProdPrec = 2.5;
  
  TVectorD v_zero(v_deltaV.GetNoElements());
  
  v_deltaV *= -1.0;
  v_deltaV += v_vdepl(getClosestIndex(v_sector,sector));
  
  TVectorD v_hig_Draw(v_deltaV);
  v_hig_Draw += v_hig;
  v_hig_Draw += vdeplProdPrec;
  
  TVectorD v_low_Draw(v_deltaV);
  v_low_Draw -= v_low;
  v_low_Draw -= vdeplProdPrec;
  
  
  std::vector<double> v_time_val;
  std::vector<double> v_time_err;
  std::vector<double> v_vdepl_val;
  std::vector<double> v_vdepl_sta;
  std::vector<double> v_vdepl_tot;
  
  // Extract for each fill the measurement and the corresponding time
  // and fill them into a STL vector
  for (int i=0; i<v_fill.GetNoElements(); ++i) {
    double time_val   = v_time(i);
    double time_err   = 0.0;
    TString s_rad("");
    TString vn_data(Form("%s/%s/%d/%d/v_volt_val%d",
                         det.Data(),lay.Data(),sector,(int)v_fill(i),val));
    if(radius>0.0) {
      s_rad = Form("%d",(int)radius);
      vn_data         += "_r"+s_rad;
    }
    
    TVectorD* vp_data = (TVectorD*)f_data->Get(vn_data.Data());

    if (!vp_data){
      Warning("plotSector","No vector found for %s",vn_data.Data());
      continue;
    }

    // Exclude fill 1616 -> used for Vdepl normalisation
    if((int)v_fill(i)==1616){
      continue;
    }
    
    TVectorD v_data   = *vp_data;
    double vdepl_val  = v_data(0);
    double vdepl_sta  = v_data(1);
    double vdepl_sys  = v_data(2);
    
    if (det.EqualTo("IT")) {
      vdepl_sys = TMath::Sqrt(vdepl_sys*vdepl_sys+STTool::voltSyst_IT*STTool::voltSyst_IT);
    }
    if (det.EqualTo("TT")) {
      vdepl_sys = TMath::Sqrt(vdepl_sys*vdepl_sys+STTool::voltSyst_TT*STTool::voltSyst_TT);
    }
    
    if (vdepl_sta>20.0) {
      Info("plotSector","Statistical uncertainty (%5.2f) for %s too high",vdepl_sta,vn_data.Data());
      continue;
    }
    v_time_val.push_back(time_val);
    v_time_err.push_back(time_err);
    v_vdepl_val.push_back(vdepl_val);
    v_vdepl_sta.push_back(vdepl_sta);
    v_vdepl_tot.push_back(TMath::Sqrt(vdepl_sta*vdepl_sta+vdepl_sys*vdepl_sys));
    
  }
  
  // Convert STL vectors into ROOT vectors
  TVectorD vr_time_val     = getROOTVector<double>(v_time_val);
  TVectorD vr_time_err     = getROOTVector<double>(v_time_err);
  TVectorD vr_vdepl_val    = getROOTVector<double>(v_vdepl_val);
  TVectorD vr_vdepl_sta    = getROOTVector<double>(v_vdepl_sta);
  TVectorD vr_vdepl_tot    = getROOTVector<double>(v_vdepl_tot);
  
  if (v_time_val.size()==0) {
    Warning("plotSector","For sector %d %s there is no data point!",sector,radius>0.0?Form("and %d mm",(int)radius):"");
    return EXIT_SUCCESS;
  }
  
  // Define graphs representing Hamburg model
  TGraph* g_hamburg = new TGraph(v_time_H,v_deltaV);
  TGraphAsymmErrors* ge_hamburg = new TGraphAsymmErrors(v_time_H,v_deltaV,v_zero,v_zero,v_low,v_hig);
  TGraph* g_hamburg_hig = new TGraph(v_time_H,v_hig_Draw);
  TGraph* g_hamburg_low = new TGraph(v_time_H,v_low_Draw);
  
  ge_hamburg->SetFillColor(kGray+2);
  g_hamburg_hig->SetLineStyle(9);
  g_hamburg_low->SetLineStyle(9);
  
  
  TGraphErrors* ge_vdepl_sta  = new TGraphErrors(vr_time_val,vr_vdepl_val,
                                                 vr_time_err,vr_vdepl_sta);
  ge_vdepl_sta->SetLineColor(kRed);
  ge_vdepl_sta->SetMarkerColor(kRed);
  ge_vdepl_sta->SetMarkerSize(1);
  
  TGraphErrors* ge_vdepl_tot  = new TGraphErrors(vr_time_val,vr_vdepl_val,
                                                 vr_time_err,vr_vdepl_tot);
  ge_vdepl_tot->SetLineColor(kRed);
  ge_vdepl_tot->SetMarkerColor(kRed);
  ge_vdepl_tot->SetMarkerSize(1);
  
  TMultiGraph* mg_frame = new TMultiGraph("mg_frame","Frame multigraph");
  
  // Plot all graphs
  TCanvas* c_frame = new TCanvas("c_frame","Frame Canvas",800,600);
  mg_frame->Draw("");
  
  mg_frame->Add(ge_hamburg,"E3");
  mg_frame->Add(g_hamburg_hig,"L");
  mg_frame->Add(g_hamburg_low,"L");
  mg_frame->Add(g_hamburg,"L");
  mg_frame->Add(ge_vdepl_sta,"P");
  mg_frame->Add(ge_vdepl_tot,"[]");
  
  mg_frame->Draw("A");
  
  // mg_frame->GetXaxis()->SetLimits(1262304000,1372636800);
  double jan2010 = 1262304000.0;
  double dec2015 = 1450103679.0;
  mg_frame->GetXaxis()->SetLimits(jan2010, dec2015);

  mg_frame->SetMinimum(100.0);
  mg_frame->SetMaximum(300.0);

  if (det.EqualTo("IT")) {
    mg_frame->SetMinimum(  0.0);
    mg_frame->SetMaximum(200.0);
  }
  mg_frame->GetXaxis()->SetTimeDisplay(1);
  mg_frame->GetXaxis()->SetTimeOffset(0.0);
  mg_frame->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%Y}");
  mg_frame->GetXaxis()->SetTitle("Date");
  mg_frame->GetYaxis()->SetTitle("#it{V}_{depl} [V]");
  mg_frame->GetXaxis()->SetLabelOffset(0.02);
  mg_frame->GetXaxis()->SetTimeOffset(0.0);
  mg_frame->GetYaxis()->SetTitleOffset(1.3);
  mg_frame->GetXaxis()->SetTitleOffset(1.4);
  
  mg_frame->Draw("A");
  
  mg_frame->GetHistogram()->Draw("axissame");
  
  if (!notSave) {
    // Save picture
    gSystem->Exec(Form("mkdir -p %s",dn_graph.Data()));
    gSystem->Exec(Form("rm %s &> /dev/null",fn_graph.Data()));
    c_frame->SaveAs(fn_graph.Data());
  }
  
  
  
  Printf("=========================================\n");

  
}


