//  
//  hamburgTool.h
//  Tools and data for Hamburg model Prediction
//
//  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"

typedef std::map<std::string,double> HMap;

namespace STHamburg {
  
  // Parameters of Hamburg model
  static std::map<std::string,double> createHamburgPara(){
    std::map<std::string,double> mp_para;
    mp_para["q"]       =     1.62e-19;
    mp_para["d"]       =   500.0e-4;
    mp_para["eps0"]    =     8.854e-14;
    mp_para["eps"]     =    11.9;
    mp_para["kB"]      =     8.617e-5;
    mp_para["Eaa"]     =     1.086;
    mp_para["Ey"]      =     1.3252;
    mp_para["ga"]      =     1.4e-2;
    mp_para["gc"]      =     1.6e-2;
    mp_para["gy"]      =     5.7e-2;
    mp_para["Nc"]      =     7.5e-2/2.29e-13;
    mp_para["c"]       =     2.29e-13;
    mp_para["ta"]      =     0.0;
    mp_para["ty"]      =     0.0;
    mp_para["k0a"]     =     2.4e13;
    mp_para["k0y"]     =     1.5e15;
    mp_para["Ny"]      =     0.0;
    return mp_para;
  }
  std::map<std::string,double> hamburgPara = createHamburgPara();
  

  // Parameter uncertainties of Hamburg model
  static std::map<std::string,double> createHamburgUnce(){
    std::map<std::string,double> mp_err;
    mp_err["q"]        =     0.0;
    mp_err["d"]        =     0.0;
    mp_err["eps0"]     =     0.0;
    mp_err["eps"]      =     0.0;
    mp_err["kB"]       =     0.0;
    mp_err["Eaa"]      =     0.03;
    mp_err["Ey"]       =     0.03;
    mp_err["ga"]       =     0.14e-2;
    mp_err["gc"]       =     0.04e-2;
    mp_err["gy"]       =     0.09e-2;
    mp_err["Nc"]       =     0.6e-2/2.29e-13;
    mp_err["c"]        =     0.0;
    mp_err["ta"]       =     0.0;
    mp_err["ty"]       =     0.0;
    mp_err["k0a"]      =     1.0e13;
    mp_err["k0y"]      =     1.1e15;
    mp_err["Ny"]       =     0.0;
    return mp_err;
  }
  std::map<std::string,double> hamburgUnce = createHamburgUnce();
  
  
  // Calculate Veff as a function of phi_n
  std::vector<double> getVeff(double phi_n,bool isIT=false){
    if (isIT) {
      hamburgPara["d"] = 410.0e-4; // So far only 410 µm sensors considered
    }
    double conv_fac    = hamburgPara["q"]*hamburgPara["d"]*hamburgPara["d"]/(2.0*hamburgPara["eps"]*hamburgPara["eps0"]);
    double Neff        = hamburgPara["Nc"]*(1-TMath::Exp(-hamburgPara["c"]*phi_n))+hamburgPara["gc"]*phi_n;
    double Neff_unc    = (TMath::Power(hamburgUnce["gc"]*phi_n,2)+
                          TMath::Power((1-TMath::Exp(-hamburgPara["c"]*phi_n))*hamburgUnce["Nc"],2)+
                          TMath::Power(TMath::Exp(-hamburgPara["c"]*phi_n)*hamburgPara["Nc"]*hamburgUnce["c"],2));
    std::vector<double> v_ret;
    v_ret.push_back(Neff*conv_fac);
    v_ret.push_back(TMath::Sqrt(Neff_unc)*conv_fac);
    
    return v_ret;
  }
  
  
};