//  
//  voltageFit.C
//  Extraction of the depletion voltage based on a sigmoid fit
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//
//  Compile with:
//  g++ -o voltageFit voltageFit.C `root-config --glibs --cflags` -lRooFit -lRooFitCore -lEG
//

//  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/incFuncROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incGraphROOT.h"
#include "../../include/incIOROOT.h"


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


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


#include "deplTool.h"


// -d detector (TT/IT)
// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2})
// -s read-out sector number
// -f fill 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
// -ss Perform a temporary save


int voltageFit(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 = voltageFit(argc,argv);
  
  
  
  Printf("End");
  if (!runBatch) {
    theApp->Run(!runBatch);
  }
  
  
  return exit;
}


//  Executable method
int voltageFit(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 = std::getenv ("DISK");
  const char* homeVar = std::getenv ("CCEHOME");
  
  int fill      = 2083;
  TString det("TT");
  TString lay("TTaU");
  int sector    = 2634;
  int val       = 3;
  bool notSave  = false;
  bool temSave  = false;
  double radius = -10.0;
  
  
  // Reading command line arguments
  for (int i=1; i<argc; i++) {
    if (strcmp(argv[i],"-d")==0 && i+1<argc) { // Detector
      det = argv[++i];
    }else if (strcmp(argv[i],"-l")==0 && i+1<argc) { // Layer
      lay = argv[++i];
    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) { // Sector
      sector = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-f")==0 && i+1<argc) { // Fill
      fill = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-v")==0 && i+1<argc) { // Number of strips
      val = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-r")==0 && i+1<argc) { // Radius
      radius = std::atof(argv[++i]);
    }else if (strcmp(argv[i],"-ns")==0) { // Do not save
      notSave = true;
    }else if (strcmp(argv[i],"-ss")==0) { // Save it in a temporary file
      temSave = true;
    }
  }
  
  
  
  
  std::vector<double> v_volt_val;
  std::vector<double> v_volt_err;
  std::vector<double> v_adc_val;
  std::vector<double> v_adc_err;
  

  
  
  
  // Vector name where the result is stored
  TString vn_result(Form("v_volt_val%d",val));
  
  
  // Directory and file name for the stored plot
  TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d/%d",
                        diskVar,det.Data(),lay.Data(),sector,fill));
  TString fn_graph_p(Form("%s/volt_val%d",
                          dn_graph.Data(),val));
  
  // Adding radius postfix
  if (radius<0.0) {
  }else{
    vn_result  +=Form("_r%d",(int)radius);
    fn_graph_p +=Form("_r%d",(int)radius);
  }
  
  fn_graph_p += ".pdf";
  
  
  Printf("=========================================");
  Info("voltageFit","Analyze: ");
  Printf("               Sector:    %d",sector);
  if (radius>0.0) {
    Printf("               Radius:    %4.2f mm",radius);
  }
  Printf("               Fill:      %d",fill);
  Printf("               Detector:  %s",det.Data());
  Printf("               Layer:     %s",lay.Data());
  
  
  
  // Directory and filename for input data
  TString dn_data(Form("%s/data/ST/Aging"
                        ,diskVar));
  TString fn_data(Form("%s/CCEScan.root",
                       dn_data.Data()));
  TString fn_dummy(Form("%s/check",dn_data.Data()));
  
  // TString dn_ratio(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db"
  //                       ,homeVar));
  TString dn_ratio(Form("%s/data/db", homeVar));
  TString fn_ratio(Form("%s/ratio.root",
                       dn_ratio.Data()));
  
  if (temSave) {
    fn_data = TString(Form("%s/temp_%s_%s_%d_%d_r%d.root",
                           dn_data.Data(),det.Data(),lay.Data(),
                           sector,fill,(int)radius));
  }

  TFile* f_data       = TFile::Open(fn_data.Data());
  
  if (!f_data) {
    Error("voltageFit","Data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TFile* f_ratio      = TFile::Open(fn_ratio.Data());
  
  if (!f_ratio) {
    Error("voltageFit","Ratio data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  
  // Load the data
  Info("voltageFit","Loading data:");
  TString dn_input(Form("%s/%s/%d/%d",
                        det.Data(),
                        lay.Data(),
                        sector,
                        fill));
  for (int i=0; i<11; i++) {
    TString vn_input(Form("%s/v_pulse%d_val%d",dn_input.Data(),i,val));
    if (radius>0.0) {
      vn_input  +=Form("_r%d",(int)radius);
    }
    TVectorD* vp_input = (TVectorD*)f_data->Get(vn_input.Data());
    if (!vp_input) {
      Warning("voltageFit","Vector %s not available",vn_input.Data());
      continue;
    }
   
    double volt = 0.0;
    
    // Extract the voltage value
    if (det.EqualTo("TT")) {
      volt = STTool::GetTTVoltage(i*6);
    }else{
      volt = STTool::GetITVoltage(i*6);
    }
    
    if (det.EqualTo("IT") && volt>250.0){
      continue;
    }
    
    
    // Get the input values
    TVectorD v_input = (*vp_input);
    
    double height_val = v_input(0);
    double height_err = v_input(1);
    
    if (det.EqualTo("IT")) {
      height_val /= 0.1306; // transform back to the amplitude
      height_err /= 0.1306; // transform back to the amplitude
      height_val *= v_input(2)*9.0/2.0*TMath::Exp(-3.0)/(20.0);
      height_err  = height_val*TMath::Sqrt(TMath::Power(height_err/height_val,2)+
                                           TMath::Power(v_input(3)/v_input(2),2)+
                                           2*v_input(8)/(height_val*v_input(2)));
    }
    
    
    Printf("        %8.2f V: %8.4f+/-%8.4f",volt,height_val,height_err);
    
    // Load the input values in to vector
    v_adc_val.push_back(height_val);
    v_adc_err.push_back(height_err);
    v_volt_val.push_back(volt);
    v_volt_err.push_back(0.0);
    
  }
  
  if (v_adc_val.size()<3) {
    Error("voltageFit","Too few (%d) data points",v_adc_val.size());
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  // Get the corresponding root vectors
  TVectorD vr_adc_val  = getROOTVector(v_adc_val);
  TVectorD vr_adc_err  = getROOTVector(v_adc_err);
  TVectorD vr_volt_val = getROOTVector(v_volt_val);
  TVectorD vr_volt_err = getROOTVector(v_volt_err);
  
  TGraphErrors* ge_volt = new TGraphErrors(vr_volt_val,vr_adc_val,vr_volt_err,vr_adc_err);
  
  TMultiGraph* mg_volt = new TMultiGraph("mg_volt","Voltage multigraph");
  
  
  // Extract the ratio values for the sampling point
  TVectorD* vp_ratio    = (TVectorD*)f_ratio->Get(Form("v_ratio_%s",lay.Data()));
  TVectorD* vp_secRatio = (TVectorD*)f_ratio->Get(Form("v_sector_%s",lay.Data()));
  
  if (!vp_ratio || !vp_secRatio) {
    Error("voltageFit","Ratio vectors not available");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TVectorD v_ratio      = *vp_ratio;
  TVectorD v_secRatio   = *vp_secRatio;
  
  TF1* fu_volt = STTool::sigmoid4;
  //TF1* fu_volt = STTool::DBfunc;
  
  double ratio_val = STTool::ITratio_val;
  double ratio_err = STTool::ITratio_err;
  
  ge_volt->Fit(fu_volt,"M0","0");
  ge_volt->Fit(fu_volt,"M0","0");
  
  double chi2ndf = fu_volt->GetChisquare()/fu_volt->GetNDF();
  Info("voltageFit","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF());
  if (chi2ndf>1.0) {
    Info("voltageFit","scaling");
    vr_adc_err*=TMath::Sqrt(chi2ndf);
  }

  

  ge_volt = new TGraphErrors(vr_volt_val,vr_adc_val,vr_volt_err,vr_adc_err);

  ge_volt->Fit(fu_volt,"M0","0");
  
  Info("voltageFit","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF());
  
  chi2ndf = fu_volt->GetChisquare()/fu_volt->GetNDF();
  if (chi2ndf<0.5) {
    Info("voltageFit","scaling");
    vr_adc_err*=TMath::Sqrt(chi2ndf);
  }
  
  ge_volt = new TGraphErrors(vr_volt_val,vr_adc_val,vr_volt_err,vr_adc_err);
  
  ge_volt->Fit(fu_volt,"M0","0");
  
  Info("voltageFit","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF());
  
  if (strcmp(det.Data(),"TT")==0) {
    ratio_val = STTool::TTratio_val;
    ratio_err = STTool::TTratio_err;
  }
  
  // Using ratio per bin
  int ratioInd = getClosestIndex<double>(v_secRatio,(double)sector);
  if (!(ratioInd<0) && (int)v_secRatio(ratioInd)==sector && true) {
    Info("voltageFit","Using specific ratio for sector %d",sector);
    ratio_val  = v_ratio(ratioInd);
  }
  
  Info("voltageFit","Ratio: %5.3f+/-%5.3f",ratio_val,ratio_err);
  
  
  ge_volt->Fit(fu_volt,"M0Q","0");
  ge_volt->Fit(fu_volt,"M0Q","0");
  
  chi2ndf = fu_volt->GetChisquare()/fu_volt->GetNDF();
  Info("voltageFit","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF());
  if (chi2ndf>1.0) {
    Info("voltageFit","scaling");
    vr_adc_err*=TMath::Sqrt(chi2ndf);
  }
  
  ge_volt = new TGraphErrors(vr_volt_val,vr_adc_val,vr_volt_err,vr_adc_err);
  
  ge_volt->Fit(fu_volt,"M0Q","0");
  
  TVirtualFitter *fit_volt = TVirtualFitter::GetFitter();
  Double_t* cov_volt       = fit_volt->GetCovarianceMatrix();
  
  int nPar = fu_volt->GetNpar();
  
  
  // Get the value for the depletion voltage
  double tarADC_val = ratio_val;
  double tarADC_sys = ratio_err;
  
  tarADC_val *= fu_volt->GetParameter("A_{0}");
  tarADC_sys *= fu_volt->GetParameter("A_{0}");
  
  double Vdepl_val  =  fu_volt->GetX(tarADC_val,10.0,400.0);
  
  
  // Get the uncertainty on the sampling point
  double tarADC_err = 0.0;
  for (Int_t j=0; j<nPar; j++) {
    for (Int_t k=0; k<nPar; k++) {
      if (strcmp(fu_volt->GetParName(j),"A_{0}")==0 || strcmp(fu_volt->GetParName(k),"A_{0}")==0) {
        continue;
      }
      tarADC_err += cov_volt[nPar*j+k]*fu_volt->GradientPar(j,&Vdepl_val)*fu_volt->GradientPar(k,&Vdepl_val);
    }
  }
  tarADC_err = TMath::Sqrt(tarADC_err);
  
  // Translate the uncertainty to the uncertainty on the depletion voltage
  double Vdepl_err_high   =  fu_volt->GetX(tarADC_val+tarADC_err,10.0,400.0);
  double Vdepl_err_low    =  fu_volt->GetX(tarADC_val-tarADC_err,10.0,400.0);
  
  double Vdepl_sys_high   =  fu_volt->GetX(tarADC_val+tarADC_sys,10.0,400.0);
  double Vdepl_sys_low    =  fu_volt->GetX(tarADC_val-tarADC_sys,10.0,400.0);
  
  
  double Vdepl_err = TMath::Abs((Vdepl_err_high-Vdepl_err_low)/2.0);
  double Vdepl_sys = TMath::Abs((Vdepl_sys_high-Vdepl_sys_low)/2.0);

  // Draw the result
  mg_volt->Add(ge_volt,"pe");
  
  TCanvas* c_vol = new TCanvas("c_vol","Voltage Canvas",800,600);
  
  double maxX = 500.0;
  if (det.EqualTo("IT")) {
    maxX = 250.0;
  }
  
  mg_volt->Draw("A");
  mg_volt->GetXaxis()->SetLimits(0.0,maxX);
  mg_volt->GetXaxis()->SetTitle("#it{V}_{bias} [V]");
  mg_volt->SetMaximum(50.0);
  mg_volt->SetMinimum(0.0);
  mg_volt->GetYaxis()->SetTitle("Charge Equivalent [ADC value]");
  
  Int_t err_n = 103;
  Double_t err_x[err_n], err_y[err_n];
  
  for (int i=0; i<err_n-2; i++) {
    err_x[i] = Vdepl_val + (i-50.0)/50.0*TMath::Sqrt(Vdepl_err*Vdepl_err+Vdepl_sys*Vdepl_sys);
    err_y[i] = fu_volt->Eval(err_x[i]);
  }
  err_x[err_n-2] = err_x[err_n-3];
  err_x[err_n-1] = err_x[0];
  err_y[err_n-2] = 0.0;
  err_y[err_n-1] = 0.0;
  
  TPolyLine *err_poly = new TPolyLine(err_n,err_x,err_y);
  err_poly->SetFillColor(kGray);
  err_poly->SetLineColor(kGray);
  err_poly->SetLineWidth(0);
  err_poly->Draw("f");
  
  for (int i=0; i<err_n-2; i++) {
    err_x[i] = Vdepl_val + (i-50.0)/50.0*Vdepl_err;
    err_y[i] = fu_volt->Eval(err_x[i]);
  }
  err_x[err_n-2] = err_x[err_n-3];
  err_x[err_n-1] = err_x[0];
  err_y[err_n-2] = 0.0;
  err_y[err_n-1] = 0.0;
  
  err_poly = new TPolyLine(err_n,err_x,err_y);
  err_poly->SetFillColor(kGray+1);
  err_poly->SetLineColor(kGray+1);
  err_poly->SetLineWidth(0);
  err_poly->Draw("f");
  
  TLine* val_line = new TLine(Vdepl_val,0.0,Vdepl_val,fu_volt->Eval(Vdepl_val));
  val_line->Draw();
  
  Printf("Drawing");
  
  fu_volt->Draw("Same");
  mg_volt->Draw("");
  
  
  mg_volt->GetHistogram()->Draw("axissame");
  
  
  Info("voltageFit","V_{depl} = (%8.4f+/-%8.4f(stat.)+/-%8.4f(syst.)) V",Vdepl_val,Vdepl_err,Vdepl_sys);
  
  
  
  // Save the result
  v_res(0) = Vdepl_val;
  v_res(1) = Vdepl_err;
  v_res(2) = Vdepl_sys;
  v_res(3) = fu_volt->GetParameter("A_{0}");
  v_res(4) = fu_volt->GetParError(fu_volt->GetParNumber("A_{0}"));
  
  while (!gSystem->Exec(Form("ls %s &> /dev/null",fn_dummy.Data()))) {
    Info("landauFit","Waiting until file is closed");
    gSystem->Sleep(1000);
  }
  gSystem->Exec(Form("touch %s &> /dev/null",fn_dummy.Data()));
  
  f_data = new TFile(fn_data.Data(),"UPDATE");

  
  
  TString dn_result = TString::Format("%s/%s/%d/%d",
                                      det.Data(),
                                      lay.Data(),
                                      sector,
                                      fill);
  
  // Get right directory
  bool isThere = f_data->GetListOfKeys()->Contains(det.Data());
  if (!isThere) {
    f_data->mkdir(det.Data());
  }
  
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s",
                                                   det.Data(),
                                                   lay.Data()));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s",
                       det.Data(),
                       lay.Data()));
  }
  
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d",
                                                   det.Data(),
                                                   lay.Data(),
                                                   sector));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s/%d",
                       det.Data(),
                       lay.Data(),
                       sector));
  }
  
  
  isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d/%d",
                                                   det.Data(),
                                                   lay.Data(),
                                                   sector,
                                                   fill));
  if (!isThere) {
    f_data->mkdir(Form("%s/%s/%d/%d",
                       det.Data(),
                       lay.Data(),
                       sector,
                       fill));
  }
  f_data->cd(Form("%s/%s/%d/%d",
                  det.Data(),
                  lay.Data(),
                  sector,
                  fill));
  if (!notSave) {
    gDirectory->Delete(Form("%s;*",
                            vn_result.Data()));
    v_res.Write(Form("%s",
                              vn_result.Data()));
    Info("voltageFit","Result written to %s",fn_data.Data());
  }
  
  f_data->Close();
  gSystem->Exec(Form("rm -f %s &> /dev/null",fn_dummy.Data()));
  if (!notSave) {
    // Save picture
    gSystem->Exec(Form("mkdir -p %s",dn_graph.Data()));
    gSystem->Exec(Form("rm %s &> /dev/null",fn_graph_p.Data()));
    c_vol->SaveAs(fn_graph_p.Data());
  }
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}