//  
//  voltageFit.C
//  Extraction of the mean ratio between the ADC value corresponding to the depletion voltage
//  and the saturation value based on a sigmoid 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 "../../include/basicROOTTools.h"


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


#include "deplTool.h"





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      = 1944;
  TString det("TT");
  TString lay("TTaU");
  int val       = 3;
  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],"-l")==0 && i+1<argc) {
      lay = argv[++i];
    }else if (strcmp(argv[i],"-v")==0 && i+1<argc) {
      val = std::atoi(argv[++i]);
    }else if (strcmp(argv[i],"-ns")==0) {
      notSave = true;
    }
  }
  
  
  
  
  
  

  
  
  

  
    
  
  Printf("=========================================");
  Info("ratioFit","Analyze: ");
  Printf("               Fill:      %d",fill);
  Printf("               Detector:  %s",det.Data());
  Printf("               Layer:     %s",lay.Data());
  
  
  
  TString dn_vdepl(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_vdepl(Form("%s/vdepl.root",
                        dn_vdepl.Data()));
  TString dn_fluka(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/fluka",homeVar));
  TString fn_fluka(Form("%s/vdepl.root",
                        dn_fluka.Data()));
  
  TFile* f_vdepl    = TFile::Open(fn_vdepl.Data());
  
  
  if (!f_vdepl) {
    Error("ratioFit","File containing V_depl does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TVectorD v_sector       = *(TVectorD*)f_vdepl->Get(Form("v_sector_%s",
                                                          lay.Data()));
  TVectorD v_vdepl        = *(TVectorD*)f_vdepl->Get(Form("v_vdepl_%s",
                                                          lay.Data()));
  
  f_vdepl->Close();
  
  Double_t a_tt_excl[18] = {0};
    /*{ 2625,2626,2627,2628,2629,2630,
      2631,2632,2633,2634,2635,2636,
      2637,2638,2639,2640,2641,2642 };*/
      
  
  TVectorD v_tt_excl(sizeof(a_tt_excl)/sizeof(a_tt_excl[0]),a_tt_excl);
  
  
  
  
  TString dn_data(Form("%s/data/ST/Aging"
                        ,diskVar));
  TString fn_data(Form("%s/CCEScan.root",
                       dn_data.Data()));

  TFile* f_data       = TFile::Open(fn_data.Data());
  
  if (!f_data) {
    Error("ratioFit","Data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TH1F* h_ratio = new TH1F("h_ratio","ratio histo",15,0.7,1.0);
  h_ratio->SetXTitle("ratio");
  h_ratio->SetYTitle(Form("fraction / (%2.2f)",0.02));
  h_ratio->Sumw2();
  
  std::vector<double> v_ratio;
  std::vector<double> v_secRatio;
  
  for (int is = 0; is<v_sector.GetNoElements(); is++) {
    int sector     = (int)v_sector(is);
    double vdepl   = v_vdepl(is);
    
    std::vector<double> v_volt_val;
    std::vector<double> v_volt_err;
    std::vector<double> v_adc_val;
    std::vector<double> v_adc_err;
    
    TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d/ratio",
                          diskVar,det.Data(),lay.Data(),sector));
    TString fn_graph_p(Form("%s/ratio_volt_val%d",
                            dn_graph.Data(),val));
    
    fn_graph_p += ".pdf";

    
    if (strcmp(det.Data(),"TT")==0) {
      bool skip = false;
      for (int ie=0; ie<v_tt_excl.GetNoElements(); ie++) {
        if ((int)v_tt_excl(ie)==sector) {
          Info("ratioFit","Excluded sector %d skipped!",sector);
          skip = true;
          break;
        }
      }
      if (skip) {
        continue;
      }
    }
    
    Info("ratioFit","Loading data for sector %d:",sector);
    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));
      TVectorD* vp_input = (TVectorD*)f_data->Get(vn_input.Data());
      if (!vp_input) {
        Warning("ratioFit","Vector %s not available",vn_input.Data());
        continue;
      }
      
      double volt = 0.0;
      
      if (det.EqualTo("TT")) {
        volt = STTool::GetTTVoltage(i*6);
      }else{
        volt = STTool::GetITVoltage(i*6);
      }
      
      if (det.EqualTo("IT") && volt>250.0) {
        continue;
      }
      
      TVectorD v_input = (*vp_input);
      
      double height_val = v_input(0);
      double height_err = v_input(1);
      
      //    height_val *= v_input(2)*v_input(2)*9.0/2.0*TMath::Exp(-3.0)/(250.0);
      //
      //    height_err  = height_val*TMath::Sqrt(TMath::Power(height_err/height_val,2)+
      //					 TMath::Power(2*v_input(3)/v_input(2),2)+
      //					 2*v_input(8)/(height_val*v_input(2)));
      
      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)));
      }

      
      //    height_val *= v_input(2)*v_input(2)*9.0/2.0*TMath::Exp(-3.0)/(250.0);
      //
      //    height_err  = height_val*TMath::Sqrt(TMath::Power(height_err/height_val,2)+
      //					 TMath::Power(2*v_input(3)/v_input(2),2)+
      //					 2*v_input(8)/(height_val*v_input(2)));
      
      
//      if (v_input(2)>20.0) {
//        height_val = v_input(0);
//        height_err = v_input(1);
//      }
      
      Printf("        %8.2f V: %8.4f+/-%8.4f",volt,height_val,height_err);
      
      
      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("ratioFit","Too few (%d) data points",v_adc_val.size());
      continue;
    }
    
    
    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");
    
    
    
    
    TF1* fu_volt = STTool::sigmoid4;
    
    
    
    ge_volt->Fit(fu_volt,"M0Q","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,"M0Q","0");
    
    TVirtualFitter *fit_volt = TVirtualFitter::GetFitter();
    Double_t* cov_volt       = fit_volt->GetCovarianceMatrix();
    
    int nPar = fu_volt->GetNpar();
    
    double adc_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;
        }
        adc_err += cov_volt[nPar*j+k]*fu_volt->GradientPar(j,&vdepl)*fu_volt->GradientPar(k,&vdepl);
      }
    }
    adc_err = TMath::Sqrt(adc_err);
    
    mg_volt->Add(ge_volt,"pe");
    
    double frac     = fu_volt->Eval(vdepl)/fu_volt->GetParameter("A_{0}");
    double frac_err = adc_err/fu_volt->GetParameter("A_{0}");
    
    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("Max. signal height [ADC value]");
    
    TLine* val_line = new TLine(vdepl,0.0,vdepl,fu_volt->Eval(vdepl));
    val_line->Draw();
    val_line = new TLine(0.0,fu_volt->Eval(vdepl),vdepl,fu_volt->Eval(vdepl));
    val_line->Draw();
    
    fu_volt->Draw("Same");
    mg_volt->Draw("");
    
    
    mg_volt->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_p.Data()));
      c_vol->SaveAs(fn_graph_p.Data());
    }
    
    if (adc_err>0.2 && !TMath::IsNaN(adc_err)) {
      h_ratio->Fill(frac,1.0/(adc_err));
    }

    Info("ratioFit","ratio: %9.6f",frac);
    
    v_ratio.push_back(frac);
    v_secRatio.push_back(sector);
    
    
    delete c_vol;
  }
  
  h_ratio->Scale(1.0/h_ratio->Integral());
  h_ratio->SetMarkerColor(kRed);
  h_ratio->SetLineColor(kRed);
  
  h_ratio->SetLabelOffset(0.02,"X");
  TCanvas* c_ratio = new TCanvas("c_ratio","Ratio plot",800,600);
  h_ratio->Draw("");
  
  Info("ratioFit","Ratio for %s: %4.3f+/-%4.3f",
       lay.Data(),h_ratio->GetMean(),h_ratio->GetRMS());
  
  TString dn_output(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db",homeVar));
  TString fn_output(Form("%s/ratio.root",
                         dn_output.Data()));
  
  TVectorD vr_ratio       = getROOTVector<double>(v_ratio);
  TVectorD vr_secRatio    = getROOTVector<double>(v_secRatio);
  
  TFile* f_output = new TFile(fn_output.Data(),"UPDATE");
  
  TString vn_ratio     = TString::Format("v_ratio_%s",lay.Data());
  TString vn_secRatio  = TString::Format("v_sector_%s",lay.Data());
  
  f_output->Delete(Form("%s;*",vn_ratio.Data()));
  f_output->Delete(Form("%s;*",vn_secRatio.Data()));
  
  vr_ratio.Write(vn_ratio.Data());
  vr_secRatio.Write(vn_secRatio.Data());
  
  f_output->Close();
  
  Info("ratioFit","Data written to %s",fn_output.Data());
  
  
  Printf("=========================================\n");

  return EXIT_SUCCESS;
}