//  
//  Hamburg.C
//  Script to calculate development
//
//  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/incDrawROOT.h"
#include "../../include/incGraphROOT.h"
#include "../../include/incHistROOT.h"
#include "../../include/incIOROOT.h"
#include "../../include/incRooFit.h"
#include "../../include/incStatROOT.h"


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


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


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


// -d detector (TT/IT)
// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2})
// -s read-out sector number
// -r consider only hit in a distance closer than r mm from the beam axis


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

double rungeKutta4(HMap, double, double, double, double,double,bool);
double GetNc(HMap, double);
double GetNa(HMap, double, double, double, double, double);
double GetNy(HMap, double, double, double, double, double);
double conv_Neff_DeltaV(HMap);

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


//  Executable method
int Hamburg(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");
  // const char* homeVar = "/afs/cern.ch/user/e/egraveri/cmtuser/STMonitoring/STAging";
  

  TString det("TT");
  TString lay("TTaU");
  double radius = -1000.0;
  int sector = 2634;
 
  
  
  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],"-r")==0 && i+1<argc) {
      radius = std::atof(argv[++i]);
    }else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
      sector = std::atoi(argv[++i]);
    }
  }
  
  
  TString vn_output = TString::Format("v_Hamburg_%s_%d",
                                      lay.Data(),sector);
  
  if (radius>0.0) {
    vn_output += Form("_r%d",(int)radius);
  }
  
  // Directory and file name for the output
  TString dn_out(Form("%s/data/ST/Aging",
                      diskVar));
  TString fn_out(Form("%s/Hamburg.root",
                      dn_out.Data()));
  
  // Directory and file name for the temperature
  TString dn_temp(Form("%s/data/temperature",homeVar));
  TString fn_temp(Form("%s/temperature.root",dn_temp.Data()));
  
  // Directory and file name for the luminosity
  TString dn_lumi(Form("%s/data/lumi",homeVar));
  //TString dn_lumi(Form("%s/data/db",homeVar));
  TString fn_lumi(Form("%s/lumi.root",dn_lumi.Data()));
  
  // // Directory and file name for the fluence
  TString dn_sector(Form("%s/data/db",homeVar));
  TString fn_sector(Form("%s/fluence.root",dn_sector.Data()));
  
  TFile* f_temp   = TFile::Open(fn_temp.Data());
  TFile* f_lumi   = TFile::Open(fn_lumi.Data());
  TFile* f_sector = TFile::Open(fn_sector.Data());
  
  
  if (!f_temp || !f_lumi || !f_sector) {
    Error("Hamburg","Some files are not available");
    return EXIT_FAILURE;
  }

  // Extract vectors to get beam periods
  TVectorF* vp_fillTime  = (TVectorF*)f_lumi->Get("v_fillTime_tot");
  TVectorF* vp_beamEner  = (TVectorF*)f_lumi->Get("v_beamEner_tot");
  TVectorF* vp_stopTime  = (TVectorF*)f_lumi->Get("v_stopTime_tot");
  TVectorF* vp_peakLumi  = (TVectorF*)f_lumi->Get("v_peakLumi_tot");
  
  if (!vp_fillTime || !vp_beamEner || !vp_stopTime || !vp_peakLumi) {
    Error("Hamburg","Luminosity vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_fillTime    = *vp_fillTime;
  TVectorD v_beamEner    = *vp_beamEner;
  TVectorD v_stopTime    = *vp_stopTime;
  TVectorD v_peakLumi    = *vp_peakLumi;
  
  
  // Extract temperature measurements
  TVectorD* vp_temp1     = (TVectorD*)f_temp->Get("v_TT_CT6_temp");
  TVectorD* vp_temp2     = (TVectorD*)f_temp->Get("v_TT_AT6_temp");
  TVectorD* vp_time1     = (TVectorD*)f_temp->Get("v_TT_CT6_time");
  TVectorD* vp_time2     = (TVectorD*)f_temp->Get("v_TT_AT6_time");
  
  if (det.EqualTo("IT")) {
    vp_temp1     = (TVectorD*)f_temp->Get("v_IT_T3_temp");
    vp_temp2     = (TVectorD*)f_temp->Get("v_IT_B3_temp");
    vp_time1     = (TVectorD*)f_temp->Get("v_IT_T3_time");
    vp_time2     = (TVectorD*)f_temp->Get("v_IT_B3_time");
  }
  
  if (!vp_temp1 || !vp_temp2 || !vp_time1 || !vp_time2) {
    Error("Hamburg","Temperature vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_temp1       = *vp_temp1;
  TVectorD v_temp2       = *vp_temp2;
  TVectorD v_time1       = *vp_time1;
  TVectorD v_time2       = *vp_time2;

  double v_time1_max = v_time1.Max();
  double v_time2_max = v_time2.Max();
  double v_fillTime_max = v_fillTime.Max();
  Info("Hamburg", "temp1, temp2 and lumi vectors finish at %f, %f, %f", v_time1_max, v_time2_max, v_fillTime_max);
  
  // Extract information about the flux in the corresponding layer
  TString vn_sector("v_sector_");
  TString vn_flux7_v("v_flux7_v_");
  TString vn_flux7_e("v_flux7_e_");
  TString vn_flux8_v("v_flux8_v_");
  TString vn_flux8_e("v_flux8_e_");
  TString vn_flux13_v("v_flux13_v_");
  TString vn_flux13_e("v_flux13_e_");
  
  vn_sector   += lay.Data();
  vn_flux7_v  += lay.Data();
  vn_flux7_e  += lay.Data();
  vn_flux8_v  += lay.Data();
  vn_flux8_e  += lay.Data();
  vn_flux13_v  += lay.Data();
  vn_flux13_e  += lay.Data();
  
  // Add radius postfix
  if (det.EqualTo("IT")) {
    
  }else{
    if (radius>0.0) {
      vn_sector   += Form("_r%d",(int)radius);
      vn_flux7_v  += Form("_r%d",(int)radius);
      vn_flux7_e  += Form("_r%d",(int)radius);
      vn_flux8_v  += Form("_r%d",(int)radius);
      vn_flux8_e  += Form("_r%d",(int)radius);
      vn_flux13_v  += Form("_r%d",(int)radius);
      vn_flux13_e  += Form("_r%d",(int)radius);
    }
  }
  
  
  
  TVectorD* vp_sector   = (TVectorD*)f_sector->Get(vn_sector.Data());
  TVectorD* vp_flux7_v  = (TVectorD*)f_sector->Get(vn_flux7_v.Data());
  TVectorD* vp_flux7_e  = (TVectorD*)f_sector->Get(vn_flux7_e.Data());
  TVectorD* vp_flux8_v  = (TVectorD*)f_sector->Get(vn_flux8_v.Data());
  TVectorD* vp_flux8_e  = (TVectorD*)f_sector->Get(vn_flux8_e.Data());
  TVectorD* vp_flux13_v  = (TVectorD*)f_sector->Get(vn_flux13_v.Data());
  TVectorD* vp_flux13_e  = (TVectorD*)f_sector->Get(vn_flux13_e.Data());

  if (!vp_sector || !vp_flux7_v || !vp_flux7_e || !vp_flux8_v || !vp_flux8_e || !vp_flux13_v || !vp_flux13_e) {
    Error("Hamburg","Flux vectors are not available");
    return EXIT_FAILURE;
  }
  
  TVectorD v_sector   = *vp_sector;
  TVectorD v_flux7_v  = *vp_flux7_v;
  TVectorD v_flux7_e  = *vp_flux7_e;
  TVectorD v_flux8_v  = *vp_flux8_v;
  TVectorD v_flux8_e  = *vp_flux8_e;
  TVectorD v_flux13_v  = *vp_flux13_v;
  TVectorD v_flux13_e  = *vp_flux13_e;
  
  // Define number of calculated pathes
  int nErr = 1000;
  
  // Parameters for Hamburg model
  HMap para_base = STHamburg::hamburgPara;
  HMap unce_base = STHamburg::hamburgUnce;
  
  if (det.EqualTo("IT")) {
    para_base["d"] = 410.0e-4; // So far only 410 µm sensors considered
  }
  
  std::vector<HMap> v_para;
  std::vector<double> v_tempUnce;
  std::vector<double> v_fluxUnce7;
  std::vector<double> v_fluxUnce8;
  std::vector<double> v_fluxUnce13;
  
  std::vector<double> v_deltaV;
  std::vector<TVectorD> v_N;
  
  std::vector<double> v_time;
  std::vector<double> v_val;
  std::vector<double> v_low;
  std::vector<double> v_hig;
  
  std::vector<double> v_flux;
  
  int i_flux = getClosestIndex(v_sector,(double)sector);
  
  double flux7_v = v_flux7_v(i_flux);
  double flux7_e = v_flux7_e(i_flux);
  double flux8_v = v_flux8_v(i_flux);
  double flux8_e = v_flux8_e(i_flux);
  double flux13_v = v_flux13_v(i_flux);
  double flux13_e = v_flux13_e(i_flux);
  
  // Temperature uncertainty of 2 degree
  double temp_e  = 2.0;
  
  Info("Hamburg","Filling parameter vectors");
  for (int i=0; i<nErr; i++) {
    HMap temp_para = para_base;
    // Define parameter for each path
    if (i==0) {
      v_para.push_back(para_base);
      v_fluxUnce7.push_back(0.0);
      v_fluxUnce8.push_back(0.0);
      v_fluxUnce13.push_back(0.0);
      v_tempUnce.push_back(0.0);
    }else{
      // Make gaussian variation of parameters
      std::map<std::string, double>::iterator it = temp_para.begin();
      for(;it!=temp_para.end(); ++it){
        it->second += gRandom->Gaus()*unce_base[it->first];
      }
      v_para.push_back(temp_para);
      v_fluxUnce7.push_back(flux7_e*gRandom->Gaus());
      v_fluxUnce8.push_back(flux8_e*gRandom->Gaus());
      v_fluxUnce13.push_back(flux13_e*gRandom->Gaus());
      v_tempUnce.push_back(temp_e*(2.0*gRandom->Uniform()-1.0));
    }
    v_deltaV.push_back(0.0);
    v_flux.push_back(0.0);
    v_N.push_back(TVectorD(3));
  }
  
  double flux   = 0.0;
  double dFlux  = 0.0;
  
  int i_temp1   = 0;
  int i_temp2   = 0;
  
  int i_lumi    = 0;
  
  
  
  double tStart = 1262325600;  // 1-1-2010
  // double tStop  = 1372654800;  // 1-7-2013
  // double tStop = 1448928000; // 1-12-2015
  double tStop = 1577750400; // 31-12-2019

  std::vector<double> endTimes;
  endTimes.push_back(tStop);
  endTimes.push_back(v_time1_max);
  endTimes.push_back(v_time2_max);
  endTimes.push_back(v_fillTime_max);
  TVectorD v_endTimes(endTimes.size(), &endTimes[0]);
  tStop = v_endTimes.Min();

  double tRun   = tStart;
  double tFillEnd = tStart;
  
  double tStep  = 60000.0;
  
  time_t raw_tStart = (long)tStart;
  time_t raw_tStop  = (long)tStop;
  
  Info("Hamburg","Starting on:  %20s",ctime(&raw_tStart));
  Info("Hamburg","Stopping on:  %20s",ctime(&raw_tStop ));
  Info("Hamburg","Steps are %6.2f sec",tStep);
  Info("Hamburg","Start running");
  
  int print_interval = 0;
  int print_freq = 200;
  
  // Run between the start and the stop time and check if the index for temperature and luminosity needs to be increased
  while (tRun<tStop) {
    double temp   = 273.15;
    
    while (i_temp1<v_time1.GetNoElements() && tRun>v_time1(i_temp1)) {
      ++i_temp1;
    }
    while (i_temp2<v_time2.GetNoElements() && tRun>v_time2(i_temp2)) {
      ++i_temp2;
    }
    temp = (v_temp1(i_temp1)+v_temp2(i_temp2))/2.0+273.15;
    
    if(i_lumi<v_fillTime.GetNoElements() && tRun>v_fillTime(i_lumi)) {
      tFillEnd = v_stopTime(i_lumi);
      
    }
    
    while (i_lumi<v_fillTime.GetNoElements() && tRun>v_stopTime(i_lumi)){
      ++i_lumi;
    }
    
    // Calculate for each path the evolution of the depletion voltage
    for (int i=0; i<v_deltaV.size(); i++) {
      
      temp += v_tempUnce[i];
      dFlux = 0.0;
      if (i_lumi<v_fillTime.GetNoElements()) {
        if (tRun<v_fillTime(i_lumi)) {
          dFlux = 0.0;
        }else{
          if (v_beamEner(i_lumi)==7.0) {
            dFlux = (flux7_v+v_fluxUnce7[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
          }else if(v_beamEner(i_lumi)==8.0) {
            dFlux = (flux8_v+v_fluxUnce8[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
          }else if(v_beamEner(i_lumi)==13.0) {
            dFlux = (flux13_v+v_fluxUnce13[i])*v_peakLumi(i_lumi)*STTool::FLUKAConvFac;
          }else{
            dFlux = 0.0;
          }
        }
      }
      
      
      // Calculate change in the effective doping concentration
      v_N[i](0) = GetNc(v_para[i],v_flux[i]);
      v_N[i](1) = GetNa(v_para[i],v_N[i](1),v_flux[i],dFlux,tStep,temp);
      v_N[i](2) = GetNy(v_para[i],v_N[i](2),v_flux[i],dFlux,tStep,temp);
      
      // Convert it to a change in the depletion voltage
      v_deltaV[i] = v_N[i].Sum()*conv_Neff_DeltaV(v_para[i]);

      v_flux[i] += dFlux*tStep;
      
      // Reset temperature
      temp -= v_tempUnce[i];
      
    }
    
    v_time.push_back(tRun);
    v_val.push_back(v_deltaV[0]);
    
    std::vector<double> v_temp = v_deltaV;
    
    std::sort(v_temp.begin(),v_temp.end());
    
    
    // Extract 68% confidence interval
    v_low.push_back(TMath::Max(0.0,v_deltaV[0]-v_temp[(int)(0.158*v_temp.size())]));
    v_hig.push_back(TMath::Max(0.0,v_temp[(int)(0.842*v_temp.size())]-v_deltaV[0]));
    
    if (print_interval%print_freq == 0) Printf("Delta V = (%6.2f+%6.2f/-%6.2f) V",
                                       v_deltaV[0],
                                       TMath::Max(0.0,v_temp[(int)(0.842*v_temp.size())]-v_deltaV[0]),
                                       TMath::Max(0.0,v_deltaV[0]-v_temp[(int)(0.158*v_temp.size())]));

    tRun+=tStep;
    time_t raw_tRun = (long)tRun;

    print_interval += 1;
    if (print_interval%print_freq == 0) Info("Hamburg","Date/Time:  %s",ctime(&raw_tRun));
    
  }
  
  
  // Reduce the number of data points to be plotted
  Info("Hamburg","Reduce number of data points");
  
  std::vector<double> v_val_red;
  std::vector<double> v_time_red;
  std::vector<double> v_hig_red;
  std::vector<double> v_low_red;
  
  int redFac = (int)(86400/tStep);
  // Printf("%d",v_val.size());
  for (int i=0; i<v_val.size(); i++) {
    // Printf("%d",i);
    if(i%redFac==0){
      v_val_red.push_back(v_val[i]);
      v_time_red.push_back(v_time[i]);
      v_hig_red.push_back(v_hig[i]);
      v_low_red.push_back(v_low[i]);
    }
  }
  
  TVectorD vr_val    = getROOTVector<double>(v_val_red);
  TVectorD vr_time   = getROOTVector<double>(v_time_red);
  TVectorD vr_hig    = getROOTVector<double>(v_hig_red);
  TVectorD vr_low    = getROOTVector<double>(v_low_red);
  
  
  // Store result
  TFile* f_out = new TFile(fn_out.Data(),"UPDATE");
  
  f_out->Delete(Form("%s_deltaV;*",vn_output.Data()));
  f_out->Delete(Form("%s_time;*"  ,vn_output.Data()));
  f_out->Delete(Form("%s_hig;*"   ,vn_output.Data()));
  f_out->Delete(Form("%s_low;*"   ,vn_output.Data()));
  
  vr_val.Write(Form("%s_deltaV",vn_output.Data()));
  vr_time.Write(Form("%s_time",vn_output.Data()));
  vr_hig.Write(Form("%s_hig",vn_output.Data()));
  vr_low.Write(Form("%s_low",vn_output.Data()));
  
  f_out->Close();
  
  Info("Hamburg","Data written to %s",fn_out.Data());

  Printf("=========================================\n");
  return EXIT_SUCCESS;
}


// Runge-Kutta 4th order algorithm for Hamburg model (not used)
double rungeKutta4(HMap para, double flux, double dFlux, double temp, double dT,double T_noBeam,bool isBeam){
  double dNc = 0.0;
  double dNa = 0.0;
  double dNy = 0.0;
  
  double tau_a = 1.0/(para["k0a"]*TMath::Exp(-para["Eaa"]/(para["kB"]*temp)));
  double tau_y = 1.0/(para["k0y"]*TMath::Exp(-para["Ey" ]/(para["kB"]*temp)));
  
  Printf("tau_a: %f,tau_y: %f",tau_a,tau_y);
  Printf("flux: %f, dFlux: %f, t: %f",flux,dFlux,T_noBeam);
  

    dNc += (para["gc"]+para["Nc"]*para["c"]*TMath::Exp(-para["c"]*flux))*dFlux*dT/2.0;
    dNa += (flux*para["ga"]*-1.0/tau_a+para["ga"]*dFlux)*TMath::Exp(-T_noBeam/tau_a)*dT/2.0;
    dNy += (dFlux*para["gy"]*(1.0-1.0/(1.0+T_noBeam/tau_y))+
            flux*para["gy"]/TMath::Power(1.0+T_noBeam/tau_y,2)/tau_y)*dT/2.0;
  
  flux += dFlux*dT/2.0;

    dNc += (para["gc"]+para["Nc"]*para["c"]*TMath::Exp(-para["c"]*flux))*dFlux*dT/2.0;
    dNa += (flux*para["ga"]*-1.0/tau_a+para["ga"]*dFlux)*TMath::Exp(-T_noBeam/tau_a)*dT/2.0;
    dNy += (dFlux*para["gy"]*(1.0-1.0/(1.0+T_noBeam/tau_y))+
            flux*para["gy"]/TMath::Power(1.0+T_noBeam/tau_y,2)/tau_y)*dT/2.0;

  Printf("dNc: %f,dNa: %f,dNy: %f",dNc,dNa,dNy);
  return dNc+dNa+dNy;
  
  
}

// Calculate stable damage part
double GetNc(HMap para, double flux){
  return para["Nc"]*(1.0-TMath::Exp(-para["c"]*flux))+para["gc"]*flux;
}

// Calculate change in the annealing part
double GetNa(HMap para, double Na, double flux, double dFlux, double dT, double temp){
  
  double tau_a = 1.0/(para["k0a"]*TMath::Exp(-para["Eaa"]/(para["kB"]*temp)));
  
  Na += (para["ga"]*dFlux-Na/tau_a)*dT/2.0;
  flux += dFlux*dT/2.0;
  Na += (para["ga"]*dFlux-Na/tau_a)*dT/2.0;
  
  return Na;
}

// Calculate change in the reverse annealing part
double GetNy(HMap para, double Ny, double flux, double dFlux, double dT, double temp){
  
  double tau_y = 1.0/(para["k0y"]*TMath::Exp(-para["Ey" ]/(para["kB"]*temp)));
  
  Ny += (TMath::Power((flux==0.0?0.0:Ny/(flux*para["gy"]))-1,2)*flux*para["gy"]/tau_y)*dT/2.0;
  flux += dFlux*dT/2.0;
  Ny += (TMath::Power((flux==0.0?0.0:Ny/(flux*para["gy"]))-1,2)*flux*para["gy"]/tau_y)*dT/2.0;
  
  return Ny;
}

// Convert change in the effective doping concentration to change in the depletion voltage
double conv_Neff_DeltaV(HMap para){
  return para["q"]*para["d"]*para["d"]/(2.0*para["eps"]*para["eps0"]);
}