// @(#)CTDR/Electron:$Name:  $:$Id: CTDR_Electron,cpp,v 2.0 2006/11/24 16:10:45 brun Exp $
// Author: Mathieu Benoit   24/11/06

/*************************************************************************
 * . Class CTDR_Electron				                 *
 *   Written by : Mathieu Benoit                                         *
 *   m.benoit@umontreal.ca						 *
 *************************************************************************/

#include "CTDR_Electron.h"
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
using namespace std;
const double Temps_pnormal=40e-9; // shallow trapping time 
const double Temps_pprofond=4e-6; // deep trapping time
const double Temps_d=10e-9;       // shallow detrapping time
const double pi=3.14159265;	  // pi
const double constrep=3*1.602e-19/(4*pi*8.854e-12*12.0);// repulsion constant

//////////////////////////////////////////////////////////////////////////
// CTDR_Electron						        //
// This class represent a charge element of a gamma interaction in CZT. //              
// It contains all the function to simulate its behavior 		//
// from generation to collection 					//
//////////////////////////////////////////////////////////////////////////

ClassImp(CTDR_Electron)


//______________________________________________________________________________
CTDR_Electron::CTDR_Electron()
{	
// The constructor instanciate an electron at an arbitrary position. 
// The state of the new trap met by the charge is computed, with 
// the time left before reaching the trap. 
	x=0;
	y=0;
	z=0;
	Q=1;
	mobilite=0.1;
	D=25e-4;
	double R=(float)rand()/RAND_MAX;
	Etat=libre;
	E_collection=Libre;
	Temps_depiegeage=0;
	if( R <=pow((1 + Temps_pnormal/Temps_pprofond),-1)){
		Etat_prochain_piege=normal;
		R=(float)rand()/RAND_MAX;
		Temps_piegeage=-Temps_pnormal*log(1-R);
		//cout << Temps_piegeage <<endl;
		}
	else{
		Etat_prochain_piege=profond;
		R=(float)rand()/RAND_MAX;
		Temps_piegeage=-Temps_pprofond*log(1-R);
		};
};

//______________________________________________________________________________
void CTDR_Electron::SetValeurInitial(double x0,double y0, double z0, double mobilite0,double Q0)
{	
// This function sets the position, the mobility 
// and the charge of an element.
// x0,y0,z0 : position
// mobilite0 : mobility
// Q0 charge (in e units)
	x=x0;
	y=y0;
	z=z0;
	mobilite=mobilite0;
	Q=Q0;
	//cout << Q0<<endl;

};
//______________________________________________________________________________
void CTDR_Electron::SetPosition(double x0,double y0, double z0)
{
// This function set the position of a charge element.
// x0,y0,z0 : position 
	x=x0;
	y=y0;
	z=z0;
};
//______________________________________________________________________________
void CTDR_Electron::SetEtat(Etats etat_voulu,E_piege etat_piege_voulu)
{
// This function modify the state of the charge element.
// etat_voulu : trapped or free
// etat_piege_voulu : next trap is deep or shallow

	Etat=etat_voulu;
	Etat_prochain_piege=etat_piege_voulu;
};
//______________________________________________________________________________
void CTDR_Electron::Collecte(double Z)
{
// This function set this charge element as collected,
// Z : The z position where we want to put the Electron when collect 
	z=Z;
     	E_collection = collecte;
     	Etat=piege;
     	Etat_prochain_piege=profond;
};
//______________________________________________________________________________
E_collecte CTDR_Electron::GetEcollecte()
{
// This function return the state of the charge, free or collected
        if(Etat==piege && Etat_prochain_piege==profond)
           {
                return collecte;
           }
        else{
                return Libre;
                };
};
//______________________________________________________________________________
long double CTDR_Electron::GetX()
{
// Access method for position X
	return x;
};
//______________________________________________________________________________
long double CTDR_Electron::GetY()
{
// Access method for position Y
	return y;
};
//______________________________________________________________________________
long double CTDR_Electron::GetZ()
{
// Access method for position Z
	return z;
};
//______________________________________________________________________________
long double CTDR_Electron::GetQ()
{
// Access method for charge
	return Q;
};
//______________________________________________________________________________
Etats CTDR_Electron::GetEtat()
{
// Access method for state (free or trapped)
	return Etat;
};
//______________________________________________________________________________
void CTDR_Electron::Piege()
{
// This method traps the charge and determine the time it will remain trapped
	Temps_piegeage=1;
	Etat=piege;
	if (Etat_prochain_piege=profond){
		Temps_depiegeage=1000;
		}
	else {  double R=(float)rand()/RAND_MAX;
		Temps_depiegeage=-Temps_d*log(1-R);
	        };
};
//______________________________________________________________________________
void CTDR_Electron::Depiege()
{
//This method untrap a charge and determine the state of the next trap.
//The time left until meeting this trap is also computed 
	double R=(float)rand()/RAND_MAX;
	Etat=libre;
	Temps_depiegeage=0;
	if( R <=pow((1 + Temps_pnormal/Temps_pprofond),-1)){
		Etat_prochain_piege=normal;
		R=(float)rand()/RAND_MAX;
		Temps_piegeage=-Temps_pnormal*log(1-R);
		}
	else{
		Etat_prochain_piege=profond;
		R=(float)rand()/RAND_MAX;
		Temps_piegeage=-Temps_pprofond*log(1-R);
		};
};

//______________________________________________________________________________
void CTDR_Electron::Diffusion_Repulsion(bool diffusion,bool repulsion,double xrms, 
 double yrms, double zrms,int N,double dt)
{
// This method simulate the process of repulsion and diffusion for one charge
// diffusion and repulsion are tokens to determine which process to include 0=off,1=on
// xrms,yrms,zrms are the size in rms of the charge could to which the charge belong
// N is the number of charges to emulate
// dt is the time step. 
// The two process are simulated by a random walk ,using a step distributed normally, 
// with a sigma adjusted according to a modified diffusion coefficient, wich include the
// effect of repulsion , see documentation, for more details on the diffusion-repulsion formula
	if(diffusion==0 && repulsion==0){}
	else if(diffusion==1 && repulsion==0){
		double R=(float)rand()/RAND_MAX;
                double R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ x+=sqrt(2*D*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ y+=sqrt(2*D*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
		  if(R2<1 & R2>0){ z+=sqrt(2*D*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  }
	else if(diffusion==1 && repulsion==1){
	          double R=(float)rand()/RAND_MAX;
                  double R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ x+=sqrt(2*(D+mobilite*N*constrep/(15*sqrt(5.0)*xrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ y+=sqrt(2*(D+mobilite*N*constrep/(15*sqrt(5.0)*yrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
		  if(R2<1 & R2>0){ z+=sqrt(2*(D+mobilite*N*constrep/(15*sqrt(5.0)*zrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  }

	else if(diffusion==0 && repulsion==1){
		  double R=(float)rand()/RAND_MAX;
                  double R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ x+=sqrt(2*(mobilite*N*constrep/(15*sqrt(5.0)*xrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
                  if(R2<1 & R2>0){ y+=sqrt(2*(mobilite*N*constrep/(15*sqrt(5.0)*yrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  R=(float)rand()/RAND_MAX;
                  R2=(float)rand()/RAND_MAX;
		  if(R2<1 & R2>0){ z+=sqrt(2*(mobilite*N*constrep/(15*sqrt(5.0)*zrms))*dt)*sin(2*pi*R)*sqrt(-2*log(R2));};
                  }

	else{};
};


//______________________________________________________________________________
void  CTDR_Electron::Transport(double Ex, double Ey, double Ez,double dt)
{
// This function transport using Euler integration, for field (Ex,Ey,Ez), 
// considered constant over time dt. The movement equation are those 
// of charges in semi-conductors, sx= mu*E*dt;
                  double dx,dy,dz;
                  dx=-mobilite*Ex*dt; // v=mobilité * Ex   dx=vdt Ex=cte sur dx
                  dy=-mobilite*Ey*dt;
                  dz=-mobilite*Ez*dt;
                  if(dx>1e-4){dx=1e-4;};
                  if(dy>1e-4){dy=1e-4;};
                  x+=dx;
                  y+=dy;
                  z+=dz;
		  //cout << Q << ' ' << mobilite << ' ' << dt << ' ' << Q*mobilite*Ez*dt << endl;
};

//______________________________________________________________________________
void CTDR_Electron::Iterateur(double dt,int N,double xrms,double yrms,double zrms,CTDR_Field &A)
{
// This method apply one time step to the charge represented by this class.
// Trapping, transport, diffusion repulsion and collection are considered. 
// dt is the time step
// N the number of charges in the cloud
// x,y,zrms are the sizes of the charge cloud to which the charge belongs
// A is a Field object that is used to get the value of the field at the charge's
// position

	if(Etat==libre){
		if(Temps_piegeage<=0){
			Piege();
			}
		else{
			Diffusion_Repulsion(1,1,xrms,yrms,zrms,N,dt);
			if(z>0.75*A.GetLz()){
                                A.SetChamp(x,y,z);
			                    Transport(A.GetEx(),A.GetEy(),A.GetEz(),dt);
            }
            else{
              Transport(0,0,-1100.0/A.GetLz(),dt);
              };
			Temps_piegeage-=dt;
		     };
		     }
	else{
		if(Temps_depiegeage<=0){
			Depiege();
			Diffusion_Repulsion (1,1,xrms,yrms,zrms,N,fabs(Temps_depiegeage));
			if(z>0.6*A.GetLz()){
                                A.SetChamp(x,y,z);
			                    Transport(A.GetEx(),A.GetEy(),A.GetEz(),fabs(Temps_depiegeage));
            }
            else{
              Transport(0,0,-1100.0/A.GetLz(),fabs(Temps_depiegeage));
              };

			};
		Temps_depiegeage-=dt;
			};


};