//  
//  deplTool.h
//  Tools and data for depletion voltage analysis
//
//  Created by Christian Elsasser on 01.05.12.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <math.h>
#include <map>


#include "TROOT.h"
#include "TSpline.h"
#include "../../include/incBasicROOT.h"


namespace STTool {
  
  

  // Cut on tracks
  static TCut get_track_selection(TString det, int fill){
    // Selection from Christian, based on S/N
    TCut cu_track_TT_run1("(TrChi2/TrNDoF<3.0) &&"
                          "(GhostP<0.01+0.1/1.6*TrChi2/TrNDoF) &&"
                          "(GhostP<0.01-0.1/3.4*(TrChi2/TrNDoF-5.0))");
    TCut cu_track_IT_run1("(GhostP<0.1*TrChi2/TrNDoF) &&"
                          "(GhostP<0.2-0.1*(TrChi2/TrNDoF-5.0))");
    // Selection from Elena, based on S/sqrt(N) and S/sqrt(S+N)
    TCut sel_TT_from3479("((((TrChi2/TrNDoF)*TMath::Cos(0.0261799387799) "
                         "-GhostP*TMath::Sin(0.0261799387799)-1.35)^2/1.1^2 "
                         "+ (GhostP*TMath::Cos(0.0261799387799) "
                         "-(TrChi2/TrNDoF)*TMath::Sin(0.0261799387799)-0.04)^2/0.08^2) < 1.)");
    //TCut sel_TT_from2252("((((TrChi2/TrNDoF)*TMath::Cos(0.0471238898038) "
    //                     "-GhostP*TMath::Sin(0.0471238898038)-1.26)^2/1.0^2 "
    //                     "+ (GhostP*TMath::Cos(0.0471238898038) "
    //                     "-(TrChi2/TrNDoF)*TMath::Sin(0.0471238898038)-0.04)^2/0.07^2) < 1.)");
    //TCut sel_TT_from2084("((((TrChi2/TrNDoF)*TMath::Cos(0.0872664625997) "
    //                     "-GhostP*TMath::Sin(0.0872664625997)-1.15)^2/1.0^2 "
    //                     "+ (GhostP*TMath::Cos(0.0872664625997) "
    //                     "-(TrChi2/TrNDoF)*TMath::Sin(0.0872664625997)-0.05)^2/0.1^2) < 1.)");
    //TCut sel_TT_from0001("((((TrChi2/TrNDoF)*TMath::Cos(0.0715584993318) "
    //                     "-GhostP*TMath::Sin(0.0715584993318)-1.3)^2/1.4^2 "
    //                     "+ (GhostP*TMath::Cos(0.0715584993318) "
    //                     "-(TrChi2/TrNDoF)*TMath::Sin(0.0715584993318)-0.05)^2/0.11^2) < 1.)")
    if (det.EqualTo("TT")) {
        if (fill <= 3478) {
            return cu_track_TT_run1;
        }
        else {
            return sel_TT_from3479;
        }
    }
    else {
        return cu_track_IT_run1;
    }
  }
  
  // Define first run of fill
  static std::map<int,int> zeroRunOfFill(){
    std::map<int,int> runMap;
    runMap[1020]      =    69606;
    runMap[1616]      =    87147;
    runMap[1944]      =    95936;
    runMap[2083]      =   101357;
    runMap[2252]      =   104091;
    runMap[2472]      =   111267;
    runMap[2797]      =   120006;
    runMap[3108]      =   129494;
    runMap[3478]      =   135605;//TO BE UPDATED!!
    runMap[3820]      =   153583;
    runMap[3960]      =   156887;
    runMap[4518]      =   166262;
    runMap[4638]      =   168092;
    runMap[4643]      =   168257;
    runMap[4856]      =   173040;
    runMap[5162]      =   181247;
    runMap[5448]      =   185508;
    return runMap;
  }
  std::map<int,int> zeroRunMap = zeroRunOfFill();
  
  // Value of one time step used in sampling time change
  double timestep = 0.10417; // ns
  
  // Calibration step residual -> time step (i.e. step%6)
  static std::map<int,double> createTimeMap(){
    std::map<int,double> timeMap;

    timeMap[0]      =    0.0;
    timeMap[1]      =  -60.0;
    timeMap[2]      =  -30.0;
    timeMap[3]      =   30.0;
    timeMap[4]      =   60.0;
    timeMap[5]      =   90.0;

    return timeMap;
  }
  std::map<int,double> timeMap = createTimeMap();


  //Added by Michele Atzeni to take into account the changes in the timing steps after fill 2252
  // Calibration step residual -> time step (i.e. step%6)
  static std::map<int,double> createTimeMap_post(){
    std::map<int,double> timeMap_post;

    timeMap_post[0]      =    0.0;
    timeMap_post[1]      =  -90.0;
    timeMap_post[2]      =  -45.0;
    timeMap_post[3]      =   45.0;
    timeMap_post[4]      =   90.0;
    timeMap_post[5]      =   135.0;

    return timeMap_post;
  }
  std::map<int,double> timeMap_post = createTimeMap_post();

  
  // Calibration step residual -> TT voltage in volts (i.e. step/6)
  static std::map<int,double> createVoltMapTT(){
    std::map<int,double> voltMapTT; 
    voltMapTT[0]    =  400.0;
    voltMapTT[1]    =  350.0;
    voltMapTT[2]    =  300.0;
    voltMapTT[3]    =  250.0;
    voltMapTT[4]    =  225.0;
    voltMapTT[5]    =  200.0;
    voltMapTT[6]    =  175.0;
    voltMapTT[7]    =  150.0;
    voltMapTT[8]    =  125.0;
    voltMapTT[9]    =  100.0;
    voltMapTT[10]   =   60.0;
    voltMapTT[11]   =   20.0;  // Not used, could remove
    voltMapTT[12]   =   20.0;  // Not used, could remove
    return voltMapTT;
  }
  std::map<int,double> voltMapTT = createVoltMapTT();
  
  // Calibration step residual -> IT voltage in volts (i.e. step/6)
  static std::map<int,double> createVoltMapIT(){
    std::map<int,double> voltMapIT; 
    voltMapIT[0]    =  300.0;
    voltMapIT[1]    =  200.0;
    voltMapIT[2]    =  170.0;
    voltMapIT[3]    =  140.0;
    voltMapIT[4]    =  120.0;
    voltMapIT[5]    =  105.0;
    voltMapIT[6]    =   90.0;
    voltMapIT[7]    =   75.0;
    voltMapIT[8]    =   60.0;
    voltMapIT[9]    =   40.0;
    voltMapIT[10]   =   20.0;
    voltMapIT[11]   =   20.0;  // Not used, could remove
    voltMapIT[12]   =   20.0;  // Not used, could remove
    return voltMapIT;
  }
  std::map<int,double> voltMapIT = createVoltMapIT();
  
  
  // Ratio values
  
  double TTratio_val = 0.953; // with time correction
  double TTratio_err = 0.019;
  
  double ITratio_val = 0.956; // with time correction
  double ITratio_err = 0.017;
  
  
  // FLUKA conversion factor
  double FLUKAConvFac = 7.2*1e4;
  
  
  // Returns time step for a given calibration step
  //Modified to take into account the change in the time steps used
  double GetTimeStep(int step, int fill){
    if (step<0) {
      Error("GetTimeStep","Calibration step is smaller than zero, returning NaN");
      return sqrt(-1);
    }
    if( fill > 2252){
      return timeMap_post[step%6];
    }
    else{
      return timeMap[step%6];
    }
    
  }
  
  // Returns time in ns for a given calibration step
  //Modified to take into account the change in the time steps used
  double GetTime(int step, int fill){
    if (step<0) {
      Error("GetTime","Calibration step is smaller than zero, returning NaN");
      return sqrt(-1);
    }

    if( fill > 2252){
      return timeMap_post[step%6]*timestep;
    }
    else{
      return timeMap[step%6]*timestep;
    }
  }


  // Returns TT voltage in V for a given calibration step
  double GetTTVoltage(int step){
    if (step<0 || step/6>12) {
      Error("GetTTVoltage","Calibration step is out of reach, returning NaN");
      return sqrt(-1);
    }
    return voltMapTT[step/6];
  }
  
  // Returns IT voltage in V for a given calibration step
  double GetITVoltage(int step){
    if (step<0 || step/6>12) {
      Error("GetITVoltage","Calibration step is out of reach, returning NaN");
      return sqrt(-1);
    }
    return voltMapIT[step/6];
  }
  
  // Return calibration step for a given time step and TT voltage
  int GetTTCalibStep(double time, double voltage, int fill){

    std::map<int,double> timeMap_sw;
    if(fill>2252){
      timeMap_sw = timeMap_post;
      }
    else{
      timeMap_sw = timeMap;
    }

    std::map<int,double>::const_iterator timeEnd = timeMap_sw.end();
    std::map<int,double>::const_iterator voltEnd = voltMapTT.end();
    for (std::map<int,double>::const_iterator itTime = timeMap_sw.begin(); 
         itTime != timeEnd; itTime++) {
      for (std::map<int,double>::const_iterator itVolt = voltMapTT.begin(); 
          itVolt != voltEnd; itVolt++) {
        if (itTime->second == time && itVolt->second == voltage) {
          return itVolt->first*6+itTime->first;
        }
      }
    }
    Warning("GetTTCalibStep","No corresponding calibration step found, returning -1");
    return -1;
  }
  
  // Return calibration step for a given time step and IT voltage
  int GetITCalibStep(double time, double voltage, int fill){

    std::map<int,double> timeMap_sw;
    if(fill>2252){
      timeMap_sw = timeMap_post;
      }
    else{
      timeMap_sw = timeMap;
    }

    std::map<int,double>::const_iterator timeEnd = timeMap_sw.end();
    std::map<int,double>::const_iterator voltEnd = voltMapIT.end();
    for (std::map<int,double>::const_iterator itTime = timeMap_sw.begin(); 
         itTime != timeEnd; itTime++) {
      for (std::map<int,double>::const_iterator itVolt = voltMapIT.begin(); 
           itVolt != voltEnd; itVolt++) {
        if (itTime->second == time && itVolt->second == voltage) {
          return itVolt->first*6+itTime->first;
        }
      }
    }
    Warning("GetITCalibStep","No corresponding calibration step found, returning -1");
    return -1;
  }
  
  // Half Gaussian for Pulse
  TF1* createHalfGauss(){
    TF1* func = new TF1("halfGauss",
                        "[0]*"
                        "(TMath::Power((x-[1])/[2],2)/2.0-"
                        "TMath::Power((x-[1])/[2],3)/6.0)*"
                          "TMath::Exp(-(x-[1])/[2])",-30.0,30.0);
    func->SetParNames("A_{0}","t_{0}","#tau");
    func->SetParameters(42,-15.0,10.0);
    func->SetParLimits(0, 10.0,800.0);
    func->SetParLimits(1,-50.0, 10.0);
    func->SetParLimits(2,  1.0, 60.0);
    func->SetNpx(1000);
    return func;
  }
  TF1* halfGauss = createHalfGauss();
  
  
  // Level as V_d
  double VdeplLevel = 0.8;
  
  // Create sigmoid function 1 for Voltage
  Double_t func_Sig1(Double_t *x, Double_t *par)
  // par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
  {
    double per = par[3];
    double s   = TMath::Log(per/(1.0-per));
    double exp = par[2]/(1+TMath::Exp(-s*(x[0]-par[1])/(par[0]-par[1])));
    double lin = par[2]/2.0*(1+s/2.0*(x[0]-par[1])/(par[0]-par[1]));
    return TMath::Min(exp,lin);
  }
  
  TF1* createSigmoid1(){
    TF1* func = new TF1("Sigmoid1",func_Sig1,-100.0,800.0,5);
    func->SetParNames("V_{depl}","V_{0}","A_{0}","level");
    func->SetParameters(250.0,160.0,42,VdeplLevel);
    func->SetParLimits(0, 0.0, 300.0);
    func->SetParLimits(1,-20.0, 200.0);
    func->SetParLimits(2, 10.0,800.0);
    func->SetParLimits(4, 0.0,1.0);
    func->FixParameter(3,VdeplLevel);
    func->SetNpx(1000);
    return func;
  }
  TF1* sigmoid1 = createSigmoid1();
  
  // Create sigmoid function 2 for Voltage
  Double_t func_Sig2(Double_t *x, Double_t *par)
  // par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
  {
  double x0 = par[2]/par[0]+par[0]/(4*par[1]);
  double x1 = x0 - par[0]/(2*par[1]);
  if(x[0]<x1){
    return par[0]*x[0];
  }
  if(x[0]<x0){
    return par[2]-(x[0]-x0)*(x[0]-x0)*par[1];
  }
  return par[2];
  }
  
  TF1* createSigmoid2(){
    TF1* func = new TF1("Sigmoid2",func_Sig2,-100.0,800.0,3);
    func->SetParNames("C_{0}","B_{0}","A_{0}");
    func->SetParameters(3.0,0.5,42);
    func->SetParLimits(0, 0.0,  10.0);
    func->SetParLimits(1, 0.0,  20.0);
    func->SetParLimits(2,10.0, 100.0);
    func->SetNpx(1000);
    return func;
  }
  TF1* sigmoid2 = createSigmoid2();
  
  
  // Create sigmoid function 3 for Voltage
  Double_t func_Sig3(Double_t *x, Double_t *par)
  // par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
  {
  if(x[0]<par[0]){
    return -TMath::Sqrt(par[1]*(x[0]-par[0])*(x[0]-par[0])+par[2])+par[3]+TMath::Sqrt(par[2]);
  }
  return par[3];
  }
  
  TF1* createSigmoid3(){
    TF1* func = new TF1("Sigmoid2",func_Sig3,-100.0,800.0,4);
    func->SetParNames("X_{0}","C_{0}","B_{0}","A_{0}");
    func->SetParameters(250.0,0.5,2.0,42);
    func->SetParLimits(0, 0.0, 500.0);
    func->SetParLimits(1, 0.0, 200.0);
    func->SetParLimits(2, 0.0, 200.0);
    func->SetParLimits(3,10.0, 100.0);
    func->SetNpx(1000);
    return func;
  }
  TF1* sigmoid3 = createSigmoid3();
  
  
  // Create sigmoid function 4 for Voltage
  Double_t func_Sig4(Double_t *x, Double_t *par)
  // par[0] = V_d, par[1] = A_0, par[2] = per
  {
  double per = par[2];
  double s   = TMath::Log(per/(1.0-per));
  double V0  = 2.0*par[0]/(2.0+s);
  double exp = par[1]/(1+TMath::Exp(-s*(x[0]-V0)/(par[0]-V0)));
  double lin = par[1]/2.0*(1+s/2.0*(x[0]-V0)/(par[0]-V0));
  return TMath::Min(exp,lin);
  }
  
  TF1* createSigmoid4(){
    TF1* func = new TF1("Sigmoid4",func_Sig4,-100.0,800.0,3);
    func->SetParNames("V_{depl}","A_{0}","level");
    func->SetParameters(250.0,160.0,42,VdeplLevel);
    func->SetParLimits(0, 0.0, 300.0);
    func->SetParLimits(1, 10.0,800.0);
    func->FixParameter(2,VdeplLevel);
    func->SetNpx(1000);
    return func;
  }
  TF1* sigmoid4 = createSigmoid4();
  
  
  
  // Create function describing the ballistic deficit
  Double_t func_BD(Double_t *x, Double_t *par){
  // par[0] = V_d, par[1]= Q_0, par[2] = lambda_0, par[3] = lambda_1, par[4] = varepsilon_s
  // par[5] = d (fixed)
    
    double w = par[5]*TMath::Min(1.0,TMath::Sqrt(x[0]/par[0]));
    
    
    
    //Printf("V_d = %f",par[0]);
    
    //Printf("%f",TMath::Exp(-(par[5]-w)/par[2]));
  
    
    double y = par[1]/par[5]*TMath::Exp(-(par[5]-w)/par[2]);
    
    //Printf("y = %f",y);
    
    int nInt = 1000;
    double dt   = 1.0*w/nInt;
    double t    = dt/2.0;
    double integral = 0.0;
    
    
    
    for (int i=0; i<nInt; i++) {
      double eps = 1.0/par[5]*(TMath::Sqrt(x[0]*par[0])-t*par[0]/par[5]);
      double lam = par[2]+par[3]*TMath::Min(1.0,eps/par[4]);
      //Printf("w = %f, dt = %f, eps = %f",w,dt,eps);
      integral += TMath::Exp(-t/lam);
      t += dt;
    }
    
    
    return y*integral*dt;
  }
  
  
  TF1* createBD(){
    TF1* func = new TF1("BDFunc",func_BD,-100.0,800.0,6);
    func->SetParNames("V_{d}","Q_{0}","#lambda_{0}",
                      "#lambda_{1}","#varepsilon_{s}","d");
    func->SetParameters(200.0,40.0,500.0,1000.0,0.8,500.0);
    func->SetParLimits(0,20.0,  500.0);
    func->SetParLimits(1, 0.0,  200.0);
    func->SetParLimits(2,190.0,  210.0);
    func->SetParLimits(3,50.0, 3000.0);
    func->SetParLimits(4, 1.0,    1.5);
    func->FixParameter(5, 500.0);
    func->SetNpx(1000);
    return func;
  }
  
  
  TF1* DBfunc = createBD();
  
  
  // Constant function
  TF1* createConst(){
    TF1* func = new TF1("Constfunc","[0]",-100.0,800.0);
    func->SetParNames("A_{0}");
    func->SetParameter(0,40.0);
    func->SetParLimits(0,10.0,  500.0);;
    func->SetNpx(1000);
    return func;
  }
  
  
  TF1* Constfunc = createConst();
  
  
  
  // Generic get X
  double GetX(TObject* o,double y,double xmin, double xmax, int steps=1000){
    if (strcmp(o->IsA()->GetName(),"TF1")==0) {
      TF1* f = (TF1*)o;
      return f->GetX(y,xmin,xmax);
    }
    
    
    
    double x = xmin;
    double step = (xmax-xmin)/steps;
    TSpline* sp = (TSpline*)o;

    
    for (int i=0; i<steps; i++) {
      if ( (sp->Eval(x)<y && sp->Eval(x+step)>y) || (sp->Eval(x)>y && sp->Eval(x+step)<y) ) {
        if (TMath::Abs(sp->Eval(x)-y)==0.0) {
          return x;
        }else if (TMath::Abs(sp->Eval(x+step)-y)==0.0){
          return x+step;
        }
        return (x/TMath::Abs(sp->Eval(x)-y)+(x+step)/TMath::Abs(sp->Eval(x+step)-y))/(1.0/TMath::Abs(sp->Eval(x)-y)+1.0/TMath::Abs(sp->Eval(x+step)-y));
      }
      x+=step;
    }
    Warning("GetX","Has not found a suitable point, returning 0.0");
    return 0.0;
  }
  
  
  // List of sectors to be excluded
  bool IsExcluded(int sec){
    // create List
    std::vector<int> v_ex;
    // TT
    // Outermost regions -> low statistics
    v_ex.push_back(2593);
    v_ex.push_back(2594);
    v_ex.push_back(2595);
    v_ex.push_back(2596);
    v_ex.push_back(2597);
    v_ex.push_back(2598);
    v_ex.push_back(2599);
    v_ex.push_back(2600);
    v_ex.push_back(2601);
    v_ex.push_back(2602);
    v_ex.push_back(2603);
    v_ex.push_back(2604);
    v_ex.push_back(2605);
    v_ex.push_back(2606);
    v_ex.push_back(2607);
    v_ex.push_back(2608);
    v_ex.push_back(2665);
    v_ex.push_back(2666);
    v_ex.push_back(2667);
    v_ex.push_back(2668);
    v_ex.push_back(2669);
    v_ex.push_back(2670);
    v_ex.push_back(2671);
    v_ex.push_back(2672);
    v_ex.push_back(2673);
    v_ex.push_back(2674);
    v_ex.push_back(2675);
    v_ex.push_back(2676);
    v_ex.push_back(2677);
    v_ex.push_back(2678);
    v_ex.push_back(2679);
    v_ex.push_back(2680);
    
    return std::find(v_ex.begin(), v_ex.end(), sec)!=v_ex.end();
  }
  
  // Systematics on depletion voltage

  /* Elena
     This is supposed to be the one extracted from the difference
     between spline fit and linear fit.
     I cannot trace back how it is done.
     I think the uncertainty due to the uncertainty on the ratio
     is already good enough, so I set these to 0.
  */
  double voltSyst_TT = 0.0;  // 13.4; // V
  double voltSyst_IT = 0.0;  //  8.8; // V
  
  
  
  
  };