//  
//  voltageFit_spline.C
//  Extraction of the depletion voltage based on a spline fit
//
//  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 <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 "TExMap.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 ("HOME");
  
  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));
  
  //
  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));
  
  
  
  
  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_spline","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());
  
  
  
  
  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("/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging/data/db"));
  TString fn_ratio(Form("%s/ratio.root",
                       dn_ratio.Data()));
  
  TString fn_syst(Form("%s/systVolt.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_spline","Data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TFile* f_ratio      = TFile::Open(fn_ratio.Data());
  
  if (!f_ratio) {
    Error("voltageFit_spline","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_spline","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 || volt<30.0)){
      continue;
    }
    
    if (det.EqualTo("TT") && (volt<80.0)){
      continue;
    }
    
    
    // Get the input values
    TVectorD v_input = (*vp_input);
    
    double height_val = v_input(0);
    double height_err = v_input(1);
    
    if (true || 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_spline","Too few (%d) data points",v_adc_val.size());
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  // Add zero point
  v_adc_val.push_back(0.0);
  v_adc_err.push_back(0.5);
  v_volt_val.push_back(-1e-9);
  v_volt_err.push_back(0.0);
  
  // 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_spline","Ratio vectors not available");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TVectorD v_ratio      = *vp_ratio;
  TVectorD v_secRatio   = *vp_secRatio;
  
  TF1* fu_volt = STTool::Constfunc;
  
  double ratio_val = STTool::ITratio_val;
  double ratio_err = STTool::ITratio_err;
  
  if (det.EqualTo("TT")) {
    ratio_val = STTool::TTratio_val;
    ratio_err = STTool::TTratio_err;
  }
  
  // Using ratio per bin
  if (false && det.EqualTo("TT")) {
    int ratioInd = getClosestIndex<double>(v_secRatio,(double)sector);
    if (!(ratioInd<0) && (int)v_secRatio(ratioInd)==sector && true) {
      Info("voltageFit_spline","Using specific ratio for sector %d",sector);
      ratio_val  = v_ratio(ratioInd);
    }
  }
  
  
  
  
  // Define fit range
  double fitMax = 510.0;
  double fitMin = 280.0;
  
  if (det.EqualTo("IT")) {
    fitMax = 280.0;
    fitMin = 160.0;
  }
  
  
  ge_volt->Fit(fu_volt,"M0","0",fitMin,fitMax);
  ge_volt->Fit(fu_volt,"M0","0",fitMin,fitMax);
  
  double chi2ndf = fu_volt->GetChisquare()/fu_volt->GetNDF();
  
  // Adding 0/0 and saturation point to the data set
  v_adc_val.insert(v_adc_val.begin(),fu_volt->GetParameter("A_{0}"));
  v_adc_err.insert(v_adc_err.begin(),fu_volt->GetParError(fu_volt->GetParNumber("A_{0}")));
  v_volt_val.insert(v_volt_val.begin(),fitMax);
  v_volt_err.insert(v_volt_err.begin(),0.0);
  
  
  // Reverse vector to be compatible with TSpline
  std::reverse(v_adc_val.begin(), v_adc_val.end());
  std::reverse(v_adc_err.begin(), v_adc_err.end());
  std::reverse(v_volt_val.begin(), v_volt_val.end());
  std::reverse(v_volt_err.begin(), v_volt_err.end());
  
  // Get the new corresponding root vectors
  TVectorD vr_adc_val_corr  = getROOTVector(v_adc_val);
  TVectorD vr_adc_err_corr  = getROOTVector(v_adc_err);
  TVectorD vr_volt_val_corr = getROOTVector(v_volt_val);
  TVectorD vr_volt_err_corr = getROOTVector(v_volt_err);
  
  // Scale them
  Info("voltageFit_spline","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF());
  if (chi2ndf>0.00) {
    Info("voltageFit_spline","scaling");
    vr_adc_err_corr*=TMath::Sqrt(chi2ndf);
  }
  
  // Get new graph
  ge_volt = new TGraphErrors(vr_volt_val_corr,vr_adc_val_corr,vr_volt_err_corr,vr_adc_err_corr);
  
  
  // Create Spline with zero first and second derivate
  TSpline5* s5_volt = new TSpline5("s5_volt",ge_volt,"b2 e1 e2",0.0,0.0,0.0,0.0);
  
  s5_volt->SetLineColor(kRed);
  
  
  // 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}");
  
  // Get depletion voltage
  double Vdepl_val  =  STTool::GetX(s5_volt,tarADC_val,10.0,400.0);
  
  
  // Get the uncertainty on the sampling point
  // with 1/distance^2 weighted uncertainty on the ADC values of the two closest data points
  int iFirst  = -1;
  int iSecon  = -1;
  // Excluding first and last point as they have been artifically added
  for(int i=1; i<vr_adc_err_corr.GetNoElements()-1; i++){
    if (iFirst==-1) {
      iFirst = i;
    }else{
      if (TMath::Abs(vr_volt_val_corr(i)-Vdepl_val)<TMath::Abs(vr_volt_val_corr(iFirst)-Vdepl_val)) {
        iFirst = i;
        continue; // Successful -> do not fill second
      }
    }
    if (iSecon==-1) {
      iSecon = i;
    }else{
      if (TMath::Abs(vr_volt_val_corr(i)-Vdepl_val)<TMath::Abs(vr_volt_val_corr(iSecon)-Vdepl_val)) {
        iSecon = i;
      }
    }
  }
  
  
  double tarADC_err = vr_adc_err_corr(iFirst)/TMath::Power(vr_volt_val_corr(iFirst)-Vdepl_val,2)+vr_adc_err_corr(iSecon)/TMath::Power(vr_volt_val_corr(iSecon)-Vdepl_val,2);
  tarADC_err /= (1.0/TMath::Power(vr_volt_val_corr(iFirst)-Vdepl_val,2)+1.0/TMath::Power(vr_volt_val_corr(iSecon)-Vdepl_val,2));
  
  Printf("tarADC_err: %f",tarADC_err);
  
  // Translate the uncertainty to the uncertainty on the depletion voltage
  double Vdepl_err_high   =  STTool::GetX(s5_volt,tarADC_val+tarADC_err,10.0,400.0);
  double Vdepl_err_low    =  STTool::GetX(s5_volt,tarADC_val-tarADC_err,10.0,400.0);
  
  double Vdepl_sys_high   =  STTool::GetX(s5_volt,tarADC_val+tarADC_sys,10.0,400.0);
  double Vdepl_sys_low    =  STTool::GetX(s5_volt,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);
  
  // Alternative determination of the systematic uncertainty by taking the
  // difference between the determined depletion voltage and the crossing point of the horizontal
  // top and the line extrapolated from (0,0)
  
  double Vdepl_sys1 = Vdepl_val-fu_volt->GetParameter("A_{0}")/s5_volt->Derivative(0.0);
  
  // Alternative determination of the systematic uncertainty by taking the
  // difference between the determined depletion voltage and the extracted value from the sigmoid function
  
  TF1* fu_voltAlt = STTool::sigmoid4;
  ge_volt->Fit(fu_voltAlt,"M0","0");
  double Vdepl_sys2 = Vdepl_val-fu_voltAlt->GetX(tarADC_val,10.0,400.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] = s5_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] = s5_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,s5_volt->Eval(Vdepl_val));
  val_line->DrawClone();
  val_line = new TLine(Vdepl_val,0.0,Vdepl_val,50.0);
  val_line->SetLineStyle(2);
  val_line->DrawClone();
  
  Printf("Drawing");
  
  TF1* fu_lin = new TF1("f_lin","[0]*x",-100,800);
  fu_lin->SetParameter(0,s5_volt->Derivative(0.0));
  fu_lin->SetLineStyle(2);
  //fu_lin->DrawCopy("Same");
  
  fu_volt->SetLineColor(kBlack);
  fu_volt->SetLineStyle(2);
  fu_volt->DrawCopy("Same");
  
  s5_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);
  
  
  // Load map
  TFile* f_map = new TFile(fn_syst.Data(),"UPDATE");
  TExMap* m_lin = (TExMap*)f_map->Get(Form("m_lin_%s",lay.Data())); // Deviations from linear approach
  TExMap* m_sig = (TExMap*)f_map->Get(Form("m_sig_%s",lay.Data())); // Deviations from sigmoid approach
  if(!m_lin){
    Info("voltageFit_spline","Linear map does not exist! Creating it...");
    m_lin = new TExMap(1000);
  }
  if(!m_sig){
    Info("voltageFit_spline","Sigmoid map does not exist! Creating it...");
    m_sig = new TExMap(1000);
  }
  (*m_lin)((Long64_t)(sector+10000*fill)) = *(reinterpret_cast<Long64_t*>(&Vdepl_sys1));
  (*m_sig)((Long64_t)(sector+10000*fill)) = *(reinterpret_cast<Long64_t*>(&Vdepl_sys2));
  Info("voltageFit_spline","Deviation to linear approach:  %8.4f V",Vdepl_sys1);
  Info("voltageFit_spline","Deviation to sigmoid approach: %8.4f V",Vdepl_sys2);
  f_map->Delete(Form("m_lin_%s;*",lay.Data()));
  f_map->Delete(Form("m_sig_%s",lay.Data()));
  m_lin->Write(Form("m_lin_%s",lay.Data()));
  m_sig->Write(Form("m_sig_%s",lay.Data()));
  f_map->Close();
  
  // 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("voltageFit_spline","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_spline","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;
}