#include "Rtypes.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Trap.h"
#include "fizconst.h"
/*static double ME     =   0.51099906;           // MeV
static double Clight = 299792458;   
static double mh     = ME*0.558;
static double me     = ME*1.08;             //ev/C^2   
static double Kboltz = 8.617385e-5;           // eV/K
static double SIgap  = 1.12;                     //eV
static double hbar   = 6.5821220e-22;            //MeV s
static double hbarc  = 197.327053;                //MeV fm*/

ClassImp(Trap)

 Trap::Trap(Double_t x1,Double_t x2,Double_t x3,Double_t x4,Double_t x5,Double_t x6,Double_t x7,Double_t x8)
{
  //Constructor:
  //	Double_t x1 ; Ratio cn/cp;
  //	Double_t x2 ; Charge charge>0 holes, charge<0 electrons;
  //	Double_t x3 ; Introduction rate for effective trapp at reference temperature
  // 	Double_t x4 ; Effective energy Level of the defect
  //	Double_t x5 ; Temperature
  //	Double_t x6 ; Fluence;
  //	Double_t x7 ; RefTemperature;
  // 	Double_t x8 ; cp (used for calculation of effective detrapping time;
Cross=x1;
Charge=x2;
IntroRate=x3;
EnergyLevel=x4;
Temperature=x5;
Fluence=x6;
RefTemperature=x7;
Cp=x8;
init=1;
TrapOnly=0;
}

 Trap::~Trap()
{
Clear();
init=0;
}

 Double_t Trap::ThermalVelocity()
{
  Double_t m;
  if(Charge > 0) m=mh*1e6; else m=me*1e6; 
  return TMath::Sqrt(3*Kboltz*Temperature*TMath::Power(Clight,2)/m)*1e2; //Poprava iz dne 14.4.2000 prej je bilo 1e2
}

 Double_t Trap::DenState()
{
  Double_t m;
  if(Charge > 0) m=mh; else m=me; 
  return 2*TMath::Power(m*Kboltz*1e-6*Temperature/(2*TMath::Pi()*TMath::Power(hbarc,2)),1.5)*1e39;
}

 Double_t Trap::DeTrapTime()
{
Double_t Eg=1.12,kb=8.617385e-5;
Double_t ni=1.45e10*TMath::Exp(  Eg/(2*kb)* ( 1/300. - 1/Temperature )   );
Double_t Ei=Eg/2+0.75*kb*Temperature*TMath::Log(0.517);
Double_t ksit=TMath::Exp((EnergyLevel-Ei)/(kb*Temperature));
Double_t eprob;
 if(Charge<0) eprob=ni*Cross*Cp*TMath::Sqrt(Temperature/RefTemperature)*ksit;
 if(Charge>0) eprob=ni*Cp*TMath::Sqrt(Temperature/RefTemperature)/ksit;
return (1/eprob)*1e-9;
}

 Double_t Trap::TrapTime()
{
  // Calculates trapping time from the given parameters
  // 
 Double_t ratio;
 Double_t  Eg=1.12,kb=8.617385e-5;
 Double_t  Ei=Eg/2+0.75*kb*Temperature*TMath::Log(0.517);
 Double_t  ksit=TMath::Exp((EnergyLevel-Ei)/(kb*Temperature));
 Double_t  ksit0=TMath::Exp((EnergyLevel-Ei)/(kb*RefTemperature));
 Double_t  Pt=1/(ksit*ksit*Cross+1);
 Double_t  Pt0=1/(ksit0*ksit0*Cross+1);
 if(Charge<0) ratio=(1-Pt)/(1-Pt0); 
 if(Charge>0) ratio=(Pt)/(Pt0);
 Double_t result=IntroRate*Fluence*TMath::Sqrt(Temperature/RefTemperature)*ratio;  
 return (1/result)*1e-9;
}


 void Trap::GetInfo()
{
printf("Defect characterization:n");
Dump();
printf("Trapping time %e [s]n", TrapTime());
printf("DeTrapping time %e [s]n", DeTrapTime());
}

 void Trap::Trapping(TH1F *inh, TH1F *outh)
{
Double_t binwidth=inh->GetBinWidth(1),expt,expdt;
Int_t nbininh=inh->GetNbinsX(),i,j,k,z,imax,jmax;
Int_t mul[4]={1,1,20,1};
Double_t tmax=0,dtf,tf,binwd=0,detrap=0,trap=0,detr=0;;
Double_t tr=TrapTime();
Double_t dtr=DeTrapTime();


char *hname[4]={"tx0","tx1","tx2","tx3"};
char *hdesc[4]={"x Charge Untrapped","x Charge Trapped  ","x Charge Detrapped","x Charge trapping "};

 if(Charge>0 && init) for(i=0;i<4;i++){hname[i][1]=(char)'p'; hdesc[i][0]=(char)'+';}
 if(Charge<0 && init) for(i=0;i<4;i++){hname[i][1]=(char)'n'; hdesc[i][0]=(char)'-';} 

if(init) 
  { 
  for(i=0;i<4;i++) outh[i]=TH1F(hname[i],hdesc[i],nbininh*mul[i],0,binwidth*nbininh*mul[i]);
  init=0;
  }
  else
  for(i=0;i<4;i++) outh[i].Reset(); 

for(i=1;i<=nbininh;i++) if(inh->GetBinContent(i)!=0) {tmax=inh->GetBinCenter(i); imax=i;}
jmax=outh[2].GetNbinsX();
//printf("tmax=%e imax=%dn",tmax,imax);

 for(i=1;i<=imax;i++) {
                           outh[0].SetBinContent(i,TMath::Exp(-inh->GetBinCenter(i)/tr));
      			   outh[1].SetBinContent(i,(outh[0].GetBinContent(i)*outh[1].GetBinWidth(i)/tr));
			   
                         }
if(TrapOnly==0)
   {
 binwd=outh[2].GetBinWidth(1);
 expdt=binwd/dtr;
 expt=binwd/tr;
 // expdt=1-TMath::Exp(-binwd/dtr);
 //expt=1-TMath::Exp(-binwd/tr);
 if(expdt>1) expdt=1; if(expt>1) expt=1;  
 //printf("exptr=%e expdrt=%e tr=%e dtr=%e",1-TMath::Exp(-binwd/tr),1-TMath::Exp(-binwd/dtr),expt,expdt);

 for(j=1;j<=jmax-imax;j++)
    {
      detr=outh[1].GetBinContent(1)*expdt;
      outh[1].SetBinContent(1,outh[1].GetBinContent(1)-detr);
      outh[2].AddBinContent(j+1,detr*inh->GetBinContent(1));
   for(i=2;i<=imax;i++) 
    {
      trap=detr*expt;
      //detr-=trap;
      detrap=(outh[1].GetBinContent(i)+trap)*expdt;
      outh[1].SetBinContent(i,outh[1].GetBinContent(i)+trap-detrap);
      detr+=detrap-trap;
      outh[2].AddBinContent(j+i,(detr)*inh->GetBinContent(i));
    }
   trap=0; detrap=0; detr=0;
   }
  for(i=1;i<outh[3].GetNbinsX();i++) outh[3].SetBinContent(i,inh->GetBinContent(i)*outh[0].GetBinContent(i)+outh[2].GetBinContent(i));
   } 
 else 
  for(i=1;i<outh[3].GetNbinsX();i++) outh[3].SetBinContent(i,inh->GetBinContent(i)*outh[0].GetBinContent(i));
}


 void Trap::Trapping(TH1F *inh)
{
Double_t binwidth=inh->GetBinWidth(1),expt,expdt;
Int_t nbininh=inh->GetNbinsX(),i,j,k,z,imax,jmax;
Int_t mul[4]={1,1,20,1};
Double_t tmax=0,dtf,tf,binwd=0,detrap=0,trap=0,detr=0;;
Double_t tr=TrapTime();
Double_t dtr=DeTrapTime();
TH1F *outh= new TH1F[4];

char *hname[4]={"tx0","tx1","tx2","tx3"};
char *hdesc[4]={"x Charge Untrapped","x Charge Trapped  ","x Charge Detrapped","x Charge trapping "};

 if(Charge>0 && init) for(i=0;i<4;i++){hname[i][1]=(char)'p'; hdesc[i][0]=(char)'+';}
 if(Charge<0 && init) for(i=0;i<4;i++){hname[i][1]=(char)'n'; hdesc[i][0]=(char)'-';} 

if(init) 
  { 
  for(i=0;i<4;i++) outh[i]=TH1F(hname[i],hdesc[i],nbininh*mul[i],0,binwidth*nbininh*mul[i]);
  init=0;
  }
  else
  for(i=0;i<4;i++) outh[i].Reset(); 

for(i=1;i<=nbininh;i++) if(inh->GetBinContent(i)!=0) {tmax=inh->GetBinCenter(i); imax=i;}
jmax=outh[2].GetNbinsX();
//printf("tmax=%e imax=%dn",tmax,imax);

 for(i=1;i<=imax;i++) {
                           outh[0].SetBinContent(i,TMath::Exp(-inh->GetBinCenter(i)/tr));
      			   outh[1].SetBinContent(i,(outh[0].GetBinContent(i)*outh[1].GetBinWidth(i)/tr));
			   
                         }
if(TrapOnly==0)
   {
 binwd=outh[2].GetBinWidth(1);
 expdt=binwd/dtr;
 expt=binwd/tr;
 // expdt=1-TMath::Exp(-binwd/dtr);
 //expt=1-TMath::Exp(-binwd/tr);
 if(expdt>1) expdt=1; if(expt>1) expt=1;  
 //printf("exptr=%e expdrt=%e tr=%e dtr=%e",1-TMath::Exp(-binwd/tr),1-TMath::Exp(-binwd/dtr),expt,expdt);

 for(j=1;j<=jmax-imax;j++)
    {
      detr=outh[1].GetBinContent(1)*expdt;
      outh[1].SetBinContent(1,outh[1].GetBinContent(1)-detr);
      outh[2].AddBinContent(j+1,detr*inh->GetBinContent(1));
   for(i=2;i<=imax;i++) 
    {
      trap=detr*expt;
      //detr-=trap;
      detrap=(outh[1].GetBinContent(i)+trap)*expdt;
      outh[1].SetBinContent(i,outh[1].GetBinContent(i)+trap-detrap);
      detr+=detrap-trap;
      outh[2].AddBinContent(j+i,(detr)*inh->GetBinContent(i));
    }
   trap=0; detrap=0; detr=0;
   }
  for(i=1;i<inh->GetNbinsX();i++) inh->SetBinContent(i,inh->GetBinContent(i)*outh[0].GetBinContent(i)+outh[2].GetBinContent(i));
   } 
 else 
  for(i=1;i<inh->GetNbinsX();i++) inh->SetBinContent(i,inh->GetBinContent(i)*outh[0].GetBinContent(i));
delete [] outh;
}