#include <iostream>
#include <iomanip>

void initialize();
void execute();
void holeTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) ;
void electronTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) ;
void oldSCTDigitization(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) ; 
bool electron( double x_e, double y_e, double &vx_e, double &vy_e, double &D_e) 
;
bool hole(double x_h, double y_h, double &vx_h, double &vy_h, double &D_h)
;
double induced (int istrip, double x, double y)
;
void EField( double x, double y, double &Ex, double &Ey ) ; 
double Drift_Time( double zhit) ;
double Drift_SurfaceTime(double ysurf)  ;


double mud_e(double E) ;
double mud_h(double E) ;
void init_mud_e(double T) ;
void init_mud_h(double T) ;
void initExEyArray() ;
void initPotentialValue() ;
//------parameters given exterbnally by jobOptions ------------------
   int    m_model;
   bool   m_isSCTDigiModel;
   double m_VD;
   double m_VB;
   double m_B;
   double m_T;
   double m_transportTimeStep;
   double m_transportTimeMax;
   double m_x0;
   double m_y0;
   int    m_coutLevel;
   int    m_eventNumber;
//------parameters mostly fixed but can be changed exterbnally  ------------
   double m_bulk_depth;
   double m_strip_pitch;
   double m_depletion_depth;
   double m_y_origin_min;
//-------- parameters for e, h transport --------------------------------
   double m_kB ;
   double m_e ;
   double m_vs_e;
   double m_Ec_e ;
   double m_vs_h;
   double m_Ec_h ;
   double m_diffusion;
   double m_beta_e ;
   double m_beta_h ;
   double m_theta;
   double m_driftMobility;
   double m_tanLA;
//---------- arrays of FEM analysis -----------------------------------
   double m_PotentialValue[81][115];
   double m_ExValue150[17][115];
   double m_EyValue150[17][115];
//-------histogramme ----------------------------------

void execute() {

//  storage arrays for every 0.5 nsec
  double Q_m2[50], Q_m1[50], Q_00[50], Q_p1[50], Q_p2[50] ;

    double x0 = m_x0 / 10000. ; // [micron] to [cm] 
    double y0 = m_y0 / 10000. ; // [micron] to [cm] 

//------------- one event loop ---------------------------------
//  clear charge arrays for induced current
    for (int it=0; it<50; it++) {
      Q_m2[it]=0.; Q_m1[it]=0.; Q_00[it]=0.; Q_p1[it]=0.; Q_p2[it]=0.;
    }
//  perform electron and hole transportation starting x0 and y0
    if( !m_isSCTDigiModel) {
       holeTransport    (x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2) ;
       electronTransport(x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2) ;
    } else oldSCTDigitization(x0, y0, Q_m2, Q_m1, Q_00, Q_p1, Q_p2 );

//  printout one event 
    std::cout<<"-------------------------------------------"<<std::endl;
    std::cout<<"EventNumber = "<<m_eventNumber<<", (x0,y0)=("
    <<x0*1e4<<","<<y0*1.e4<<") [microns]"<<std::endl;
    std::cout<< std::setprecision(5)<<std::setiosflags(std::ios::fixed)
    <<"   time :         Q_m2      Q_m1      Q_00      Q_p1      Q_p2"<<std::endl;
    double qm2=0., qm1=0., q00=0., qp1=0., qp2=0.;
    for (int it = 0; it< 50; it++) {
       double time = (it +1) * 0.5 ;
       std::cout <<"  "<<time<<" ns : "
       <<std::setw(10)<<Q_m2[it]<<std::setw(10)<<Q_m1[it] 
       <<std::setw(10)<<Q_00[it]<<std::setw(10)<<Q_p1[it] 
       <<std::setw(10)<<Q_p2[it]<<std::endl;
       qm2 += Q_m2[it];
       qm1 += Q_m1[it];
       q00 += Q_00[it];
       qp1 += Q_p1[it];
       qp2 += Q_p2[it];
    }
    double qt = qm2 + qm1 + q00 + qp1 + qp2;
    std::cout <<"  TOTAL SUM   : "
    <<std::setw(10)<<qm2<<std::setw(10)<<qm1 
    <<std::setw(10)<<q00<<std::setw(10)<<qp1 
    <<std::setw(10)<<qp2<<std::setw(10)<<qt<<std::endl;
//------------- end of one event loop -------------------------
   return;
}


void initialize(){
//---------------setting basic parameters---------------------------
m_model = 2 ; // model=2 for FEM field model
m_isSCTDigiModel = false ;
m_VD = 70. ; // [Volt]
m_VB = 150. ;  // [Volt]
m_B = -2.0 ;  // [Tesla]
m_T = 273.15 ;
m_transportTimeStep = 0.25  ;
m_transportTimeMax = 25.0  ;
m_coutLevel = 0 ;   // specify printout level
////------------------usually fixed-------------------------------------
m_bulk_depth = 0.0285 ; // in [cm]
m_strip_pitch = 0.0080 ;  // in [cm]
m_T = 273.15;
//--------------------------------------------------------
   std::cout<<"EfieldModel\t"<< m_model <<"\t(default = 2)"<<std::endl;
   std::cout<<"IsSCTDigiModel\t"<<m_isSCTDigiModel<<
              "\t(default = false)"<<std::endl;
   std::cout<<"DepletionVoltage_VD\t"<< m_VD  << "\t(default = 70.)"<<std::endl;
   std::cout<<"BiasVoltage_VB\t" << m_VB  << "\t(150.)"<<std::endl;
   std::cout<<"MagneticField_B\t"<< m_B  << "\t(-2.0)"<<std::endl; 
   std::cout<<"Temperature_T\t"<< m_T  << "\t(273.15)"<<std::endl;
   std::cout<<"TransportTimeStep \t"<< m_transportTimeStep << "\t(0.2)"<<
   std::endl;
   std::cout<<"TransportTimeMax\t"<< m_transportTimeMax <<"\t(30.0)"<<std::endl;
   std::cout<<"CoutLevel\t"<< m_coutLevel  << "\t(0)"<<std::endl;
   std::cout<<"BulkDepth\t"<< m_bulk_depth  << "\t(0.0285)"<<std::endl;
   std::cout<<"StripPitch\t"<< m_strip_pitch  << "\t(0.0080)"<<std::endl;
   std::cout<<"----------------------------------------------"<<std::endl;
//------------------------------------------------
//------------------------ initialize subfunctions ------------
init_mud_h(m_T) ;
init_mud_e(m_T) ;
initExEyArray();
initPotentialValue();
return;
}

//--------------------------------------------------------------
//   drift mobility for electrons
//--------------------------------------------------------------
double mud_e(double E) {
  return m_vs_e/m_Ec_e/pow(1+pow(E/m_Ec_e,m_beta_e),1./m_beta_e);
}
//--------------------------------------------------------------
//   initialize drift mobility for electrons
//--------------------------------------------------------------
void init_mud_e(double T) {
   m_vs_e=1.53E9*pow( T,-0.87);
   m_Ec_e = 1.01*pow( T, 1.55);
   m_beta_e = 2.57E-2* pow(T,0.66);
   std::cout<<"---- parameters for electron transport -----"<<std::endl;
   std::cout<<"     m_vs_e    = "<< m_vs_e << std::endl;
   std::cout<<"     m_Ec_e    = "<< m_Ec_e << std::endl;
   std::cout<<"     m_beta_e  = "<< m_beta_e << std::endl;
return ;
}

//--------------------------------------------------------------
//   drift mobility for holes
//--------------------------------------------------------------
double mud_h(double E) {
  return m_vs_h/m_Ec_h/pow(1+pow(E/m_Ec_h,m_beta_h),1./m_beta_h);
}

//--------------------------------------------------------------
//   initialize drift mobility for holes
//--------------------------------------------------------------
void init_mud_h(double T) {
  m_vs_h = 1.62E8 * pow( T, -0.52);
  m_Ec_h = 1.24 * pow( T, 1.68);
  m_beta_h = 0.46 * pow(T,0.17);
  std::cout<<"----parameters for hole transport -----"<<std::endl;
  std::cout<<"     m_vs_h    = "<< m_vs_h << std::endl;
  std::cout<<"     m_Ec_h    = "<< m_Ec_h << std::endl;
  std::cout<<"     m_beta_h  = "<< m_beta_h << std::endl;
  return ;
}
/////////////////////////////////////////////////////////////////////
double ExValue150(int& ix, int& iy), EyValue150(int& ix, int& iy);

void initExEyArray() {
  for (int ix=0 ; ix <17 ; ix++ ){
     for (int iy=0 ; iy<115 ; iy++ ) {
       m_ExValue150[ix][iy] = ExValue150(ix,iy);
       m_EyValue150[ix][iy] = EyValue150(ix,iy);
     }
   }
}
 
double  GetPotentialValue(int& ix,int& iy);

void initPotentialValue() {
  for (int ix=0 ; ix <81 ; ix++ ){
     for (int iy=0 ; iy<115 ; iy++ ) {
      m_PotentialValue[ix][iy] = GetPotentialValue(ix,iy);
     }
  }
}
//---------------------------------------------------------------------
//  holeTransport
//---------------------------------------------------------------------
void holeTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) {
// transport holes in the bulk 
// T. Kondo, 2010.9.9
// External parameters to be specified
//     m_transportTimeMax  [nsec]
//     m_transportTimeStep [nsec]     
//     m_bulk_depth        [cm]
// Induced currents are added to
//     Q_m2[50],Q_m1[50],Q_00[50],Q_p1[50],Q_p2[50] 
 
// 
   double x = x0;  // original hole position [cm]
   double y = y0;  // original hole position [cm]
   bool   isInBulk = true;
   double t_current = 0.;
   double qstrip[5];
   double vx, vy, D;

   for (int istrip = -2 ; istrip < 3 ; istrip++) 
   qstrip[istrip+2] = induced(istrip, x, y);
   while ( t_current < m_transportTimeMax ) {
     if ( !isInBulk ) break;
     if ( !hole( x, y, vx, vy, D)) break ;
     double delta_y = vy * m_transportTimeStep *1.E-9; 
     y += delta_y;
     double dt = m_transportTimeStep;
     if ( y > m_bulk_depth ) {
        isInBulk = false ;
        dt = (m_bulk_depth - (y-delta_y))/delta_y * m_transportTimeStep;
        y = m_bulk_depth;
     }
     t_current = t_current + dt;
     x += vx * dt *1.E-9;
     double diffusion = sqrt (2.* D * dt*1.E-9);
     y += diffusion * gRandom->Gaus(0.,1.);
     x += diffusion * gRandom->Gaus(0.,1.);
     if( y > m_bulk_depth) {
         y = m_bulk_depth;
         isInBulk = false;
     }

//   get induced current by subtracting induced charges
     for (int istrip = -2 ; istrip < 3 ; istrip++) {
        double qnew = induced( istrip, x, y);
        int jj = istrip + 2;
        double dq = qnew - qstrip[jj];
        qstrip[jj] = qnew ;
        if(m_coutLevel>2 && istrip==0) std::cout<<"h:t,x,y="<<t_current
        <<" "<<x*1e4<<" "<<y*1e4<<" dq[0]="<<dq<<std::endl;
        int jt = int( (t_current+0.001) / 0.50) ;
        if(jt < 50) {
        switch(istrip) {
           case -2: Q_m2[jt] += dq ; break;
           case -1: Q_m1[jt] += dq ; break;
           case  0: Q_00[jt] += dq ; break;
           case +1: Q_p1[jt] += dq ; break;
           case +2: Q_p2[jt] += dq ; break;
        } 
        }
     }
   }  // end of hole tracing 
   if (m_coutLevel>1) 
      std::cout<<"holeTransport : x,y=("<< x0*1.e4<<","<<y0*1.e4<<")->("
      <<x*1.e4<<"," <<y*1.e4 <<") t="<<t_current<<std::endl;

   return ;
}

//---------------------------------------------------------------------
//  electronTransport
//---------------------------------------------------------------------
void electronTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) {
// transport electrons in the bulk 
// T. Kondo, 2010.9.10
// External parameters to be specified
//     m_transportTimeMax  [nsec]
//     m_transportTimeStep [nsec]     
//     m_y_origin_min      [cm]  0 except under-depleted cases
// Induced currents are added to
//     Q_m2[50],Q_m1[50],Q_00[50],Q_p1[50],Q_p2[50] 
 
// 
   double x = x0;  // original electron position [cm]
   double y = y0;  // original electron position [cm]
   bool isInBulk = true;
   double t_current = 0.;
   double qstrip[5];
   double vx, vy, D;

   for (int istrip = -2 ; istrip < 3 ; istrip++) 
   qstrip[istrip+2] = -induced(istrip, x, y);

//      note  0.004[cm] is to shift from Richard's to Taka's coordinates
   while (t_current < m_transportTimeMax) {
     if ( !isInBulk ) break;
     if ( !electron( x, y, vx, vy, D)) break ;
     double delta_y = vy * m_transportTimeStep *1.E-9; 
     y += delta_y;
     double dt = m_transportTimeStep;  // [nsec]
     if ( y < m_y_origin_min ) {
        isInBulk = false ;
        dt = (m_y_origin_min - (y-delta_y))/delta_y * m_transportTimeStep;
        y = m_y_origin_min; 
     }
     t_current = t_current + dt;
     x += vx * dt *1.E-9;
     double diffusion = sqrt (2.* D * dt*1.E-9);
     y += diffusion * gRandom->Gaus(0.,1.);
     x += diffusion * gRandom->Gaus(0.,1.);
     if( y < m_y_origin_min) {
        y = m_y_origin_min;
        isInBulk = false;
     }
//   get induced current by subtracting induced charges
     for (int istrip = -2 ; istrip < 3 ; istrip++) {
        double qnew = -induced( istrip, x, y);
        int jj = istrip + 2;
        double dq = qnew - qstrip[jj];
        qstrip[jj] = qnew ;
        if(m_coutLevel>2 && istrip==0) std::cout<<"e:t,x,y="<<t_current
          <<" "<<x*1e4<<" "<<y*1e4<<" dq[0]="<<dq<<std::endl;
        int jt = int( (t_current + 0.001) / 0.50);
        if(jt< 50) {
        switch(istrip) {
           case -2: Q_m2[jt] += dq ; break;
           case -1: Q_m1[jt] += dq ; break;
           case  0: Q_00[jt] += dq ; break;
           case +1: Q_p1[jt] += dq ; break;
           case +2: Q_p2[jt] += dq ; break;
        } 
        }
     }
   }  // end of hole tracing 
   if (m_coutLevel>1) 
      std::cout<<"elecTransport : x,y=("<< x0*1.e4<<","<<y0*1.e4<<")->("
      <<x*1.e4<<"," <<y*1.e4 <<") t="<<t_current<<std::endl;

   return ;
}

//---------------------------------------------------------------------
//  old SCT Digitization Model
//---------------------------------------------------------------------
void oldSCTDigitization(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2 ) { 
// based on old SCT Digitization model used upto 15.8.0 
// T. Kondo, 2010.9.11
// External parameters to be specified
//     m_strip_pitch
//     m_y_origin_min  [cm]  0 except under-depleted cases
// Induced currents are added to
//     Q_m2[50],Q_m1[50],Q_00[50],Q_p1[50],Q_p2[50] 
//

   double deltay = m_bulk_depth - y0 ;
   double t_drift = Drift_Time(deltay) ;
   double x_hit = x0 + m_tanLA * deltay;
   double diffusion = sqrt (2.* m_diffusion * t_drift);
   x_hit += diffusion * gRandom->Gaus(0.,1.);
   int istrip = int(x_hit/m_strip_pitch+999.0)-999;
   double x_strip = ( istrip+ 0.5) * m_strip_pitch;
   double deltax = fabs((x_hit - x_strip) / m_strip_pitch);
   double t_surf =  Drift_SurfaceTime(2.*deltax);
   double t_current =  (t_drift + t_surf) * 1.e9;
   double dq = 1.0;
   if(m_coutLevel>2) 
     std::cout<<"SCT digimodel: x,y->xfinal,istrip,dx="
        <<x0*1.e4<<","<<y0*1.e4<<"->"<<x_hit*1.e4<<" "<<istrip
        <<" "<<deltax
        << " ;t_drift,t_surf->t_current="<<t_drift<<" "<<t_surf
        <<" "<<t_current<<std::endl;
//
   int jt = int( (t_current + 0.001)/0.50);
   if(jt<50) {
      switch(istrip) {
           case -2: Q_m2[jt] += dq ; break;
           case -1: Q_m1[jt] += dq ; break;
           case  0: Q_00[jt] += dq ; break;
           case +1: Q_p1[jt] += dq ; break;
           case +2: Q_p2[jt] += dq ; break;
      }
    }

   return ;
}
//////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------
//      parameters for electron transport
//---------------------------------------------------------------
bool electron( double x_e, double y_e, double &vx_e, double &vy_e, double &D_e) 
{
double kB= 1.38E-23;       // [m^2*kg/s^2/K]
double e= 1.602E-19;       // [Coulomb]
double E, Ex, Ey, mu_e, v_e, r_e, tanLA_e, secLA, cosLA, sinLA;
EField(x_e, y_e, Ex, Ey);    // [V/cm]
if( Ey > 0.) {
   double REx = -Ex; // because electron has negative charge 
   double REy = -Ey; // because electron has negative charge
   E = sqrt(Ex*Ex+Ey*Ey);
   mu_e = mud_e(E);
   v_e = mu_e * E;
   r_e = 1.13+0.0008*(m_T-273.15);
   tanLA_e = r_e * mu_e * (-m_B) * 1.E-4;  // because e has negative charge
   secLA = sqrt(1.+tanLA_e*tanLA_e);
   cosLA=1./secLA;
   sinLA = tanLA_e / secLA;
   vy_e = v_e * (REy*cosLA - REx*sinLA)/E;
   vx_e = v_e * (REx*cosLA + REy*sinLA)/E;
   D_e = kB * m_T * mu_e/ e;  
   return true;
}  else return false;
}

//---------------------------------------------------------------
//      parameters for hole transport
//---------------------------------------------------------------
bool hole(double x_h, double y_h, double &vx_h, double &vy_h, double &D_h)
{
double kB= 1.38E-23;       // [m^2*kg/s^2/K]
double e= 1.602E-19;       // [Coulomb]
double E, Ex, Ey, mu_h, v_h, r_h, tanLA_h, secLA, cosLA, sinLA;
EField( x_h, y_h, Ex, Ey);  // [V/cm]
if( Ey > 0.) {
   E = sqrt(Ex*Ex+Ey*Ey);
   mu_h = mud_h(E);
   v_h = mu_h * E;
   r_h = 0.72 - 0.0005*(m_T-273.15);
   tanLA_h = r_h * mu_h * m_B * 1.E-4;
   secLA = sqrt(1.+tanLA_h*tanLA_h);
   cosLA=1./secLA;
   sinLA = tanLA_h / secLA;
   vy_h = v_h * (Ey*cosLA - Ex*sinLA)/E;
   vx_h = v_h * (Ex*cosLA + Ey*sinLA)/E;
   D_h = kB * m_T * mu_h/ e;  
   return true;
}  else return false;
}

//-------------------------------------------------------------------
//    calculation od induced charge using Weighting (Ramo) function
//-------------------------------------------------------------------
double induced (int istrip, double x, double y)
{
// x and y are the coorlocation of charge (e or hole)
// induced chardege on the strip "istrip" situated at the height y = d
// the center of the strip (istrip=0) is x = 0.004 [cm]
   double deltax = 0.0005, deltay = 0.00025;
   
   if ( y < 0. || y > m_bulk_depth) return 0;
   double xc = m_strip_pitch * (istrip + 0.5);
   double dx = fabs( x-xc );
   int ix = int( dx / deltax );
   if ( ix > 80 ) return 0.;
   int iy = int( y  / deltay );
   double fx = (dx - ix*deltax) / deltax;
   double fy = ( y - iy*deltay) / deltay;
   int ix1 = ix + 1;
   int iy1 = iy + 1;
   double P = m_PotentialValue[ix][iy]   *(1.-fx)*(1.-fy)
            + m_PotentialValue[ix1][iy]  *fx*(1.-fy)
            + m_PotentialValue[ix][iy1]  *(1.-fx)*fy
            + m_PotentialValue[ix1][iy1] *fx*fy ;
//   cout <<"x,y,iy="<<x<<" "<<y<<" "<<iy<<" istrip,xc,dx,ix="
//   <<istrip<<" "<<xc<<" " <<dx<<" "<<ix<<" fx,fy="<<fx
//   <<" " <<fy<< ", P="<<P<<endl;
   return P;
}

//--------------------------------------------------------------
//   Electric Field Ex and Ey
//--------------------------------------------------------------
void EField( double x, double y, double &Ex, double &Ey ) { 
// m_model == 0; uniform E-field model
// m_model == 1; flat diode model
// m_model == 2; FEM solusions
// x == 0.0040 [cm] : at the center of strip
// x must be within 0 and 0.008 [cm]

   double deltay=0.00025, deltax=0.00025  ;   // [cm]  

   Ex = 0.;
   Ey = 0.;
   if( y < 0. || y > m_bulk_depth) return;
//---------- case for FEM analysis solution  -------------------------
   if(m_model==2) {       
       int iy = int (y/deltay);
       double fy = (y-iy*deltay) / deltay;
       double xhalfpitch=m_strip_pitch/2.;
       double x999 = 999.*m_strip_pitch;
       double xx = fmod (x+x999, m_strip_pitch);
       double xxx = xx;
       if( xx > xhalfpitch ) xxx = m_strip_pitch - xx;
       int ix = int(xxx/deltax) ;
       double fx = (xxx - ix*deltax) / deltax;         
//       std::cout <<"x,y,ix,iy,fx,fy,xx,xxx= "<<x<<" "<<y<<" "<<ix
//       <<" "<<iy<<" "<<fx<<" "<<fy<<" "<<xx<<" "<<xxx<<std::endl;
       int ix1=ix+1;
       int iy1=iy+1;
       double Ex00 = 0., Ex10 = 0., Ex01 = 0., Ex11 = 0.;
       double Ey00 = 0., Ey10 = 0., Ey01 = 0., Ey11 = 0.;
 //-------- pick up appropriate data file for m_VD-----------------------
       int iVB = int(m_VB);
       switch(iVB) {
       case 150: 
         Ex00 = m_ExValue150[ix][iy];  Ex10 = m_ExValue150[ix1][iy];
         Ex01 = m_ExValue150[ix][iy1]; Ex11 = m_ExValue150[ix1][iy1];
         Ey00 = m_EyValue150[ix][iy];  Ey10 = m_EyValue150[ix1][iy];
         Ey01 = m_EyValue150[ix][iy1]; Ey11 = m_EyValue150[ix1][iy1];
         break;
         } 
//---------------- end of data bank search---
       Ex = Ex00*(1.-fx)*(1.-fy) + Ex10*fx*(1.-fy)
              + Ex01*(1.-fx)*fy + Ex11*fx*fy ;
//     std::cout <<"xx, xhalfpitch = "<< xx << " "<< xhalfpitch<<std::endl;
       if(xx > xhalfpitch ) Ex = -Ex;
       Ey = Ey00*(1.-fx)*(1.-fy) + Ey10*fx*(1.-fy)
              + Ey01*(1.-fx)*fy + Ey11*fx*fy ;
       return;
   }

//---------- case for uniform electriv field ------------------------
   if( m_model ==0 ) {
       Ey = m_VB / m_depletion_depth ;   
       return;
   }
//---------- case for flat diode model ------------------------------
   if(m_model==1) {       
       if(m_VB > m_VD) 
          Ey = (m_VB+m_VD)/m_depletion_depth 
                - 2.*m_VD*(m_bulk_depth-y)/(m_bulk_depth*m_bulk_depth);
       else {
          double Emax = 2.* m_depletion_depth * m_VD /
                        (m_bulk_depth*m_bulk_depth);
          Ey = Emax*(1-(m_bulk_depth-y)/m_depletion_depth);
       }
       return;
   }

return;
}


//----------------------------------------------------------------------
// perpandicular Drift time calculation
//----------------------------------------------------------------------
double Drift_Time( double zhit) {
   double denominator = m_VD + m_VB - (2.0*zhit*m_VD/m_bulk_depth);
   double driftTime = log((m_VD + m_VB)/denominator) ;
   driftTime *= m_bulk_depth * m_bulk_depth / (2.0 * m_driftMobility * m_VD);
   return driftTime ;
}

//----------------------------------------------------------------------
// Surface Drift time calculation
//----------------------------------------------------------------------
double Drift_SurfaceTime(double ysurf)  {
   double surfaceDriftTime ;
   if(ysurf<0.5) {
      double y = ysurf/ 0.5;
      surfaceDriftTime =  5.0e-9 * y * y;
   }
   else {
      double y = (1.0 - ysurf)/ 0.5 ;
      surfaceDriftTime = 1.e-8 + ( 5.e-9 - 1.e-8) * y * y;
   }
   return surfaceDriftTime ;
}
/////////////////////////////////////////////////////////////////////