//#include "AthenaKernel/IAtRndmGenSvc.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGauss.h"
#include "SCT2/MySCTcharge_v2.h"
#include <TRandom.h>
///////////////////////////////////////----------------

MySCTcharge::MySCTcharge(const std::string& name, ISvcLocator* pSvcLocator) :
  AthAlgorithm(name, pSvcLocator),
  m_rndmSvc("AtRndmGenSvc",name),
  m_rndmEngineName("MySCTcharge")
//  Algorithm(name, pSvcLocator),m_analysisTools( "AnalysisTools", this )
//  , m_rndmEngineName(0)
{
//---------------setting basic parameters---------------------------
   declareProperty("EfieldModel", m_model = 2); // model=2 for FEM field model
   declareProperty("DepletionVoltage_VD", m_VD = 70.);   // [Volt]
   declareProperty("BiasVoltage_VB", m_VB = 150.);  // [Volt]
   declareProperty("MagneticField_B", m_B = -2.0);  // [Tesla]
   declareProperty("Temperature_T", m_T = 273.15);
   declareProperty("TransportTimeStep", m_transportTimeStep = 0.25 );
   declareProperty("TransportTimeMax", m_transportTimeMax = 25.0 );
   declareProperty("X0", m_x0 = 5.0) ;  // [micron]
   declareProperty("Y0", m_y0 = 20.0) ; // [micron]
   declareProperty("CoutLevel", m_coutLevel = 0);   // specify printout level
//------------------usually fixed-------------------------------------
   declareProperty("BulkDepth", m_bulk_depth = 0.0285); // in [cm]
   declareProperty("StripPitch", m_strip_pitch = 0.0080);  // in [cm]
}

///////////////////////////////////////////////////////////////////////////////

MySCTcharge::~MySCTcharge() {}

//////////////////////////////////////////////////////////////////////////////

StatusCode MySCTcharge::initialize() {

////////////////////////////////////////////////////////////
   std::cout<<"================ input parameters=================="<<std::endl;
   std::cout<<"EfieldModel\t"<< m_model <<"\t(default = 2)"<<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;
///////////////////////////////////////////////////////////
  m_eventNumber = 0;
  initTransport();
//+++ Get the random generator engine
  m_sc = initRandomEngine();
  if(m_sc.isFailure()) return StatusCode::FAILURE ;

  return m_sc;
}

///////////////////////////////////////////////////////////////////

StatusCode MySCTcharge::execute() {

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

    m_eventNumber++ ;
    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
       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) ;

//  printout one event 
    std::cout<<
    "=========================start of event==============================="
    <<std::endl;
    std::cout<<"EventNumber = "<<m_eventNumber<<", initial eh pair at (x0,y0)=("
    <<x0*1e4<<","<<y0*1.e4<<") [microns]"<<std::endl;
    std::cout<< std::setprecision(4)<<std::setiosflags(std::ios::fixed)<<
    "=======================induced charge================================="
    <<std::endl;
    std::cout<<
    "   time :        Q_m2     Q_m1     Q_00     Q_p1     Q_p2   :   Total"
    <<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 ;
       double sum = Q_m2[it]+Q_m1[it]+Q_00[it]+Q_p1[it]+Q_p2[it];
       std::cout <<"  "<<time<<" ns : "
       <<std::setw(9)<<Q_m2[it]<<std::setw(9)<<Q_m1[it] 
       <<std::setw(9)<<Q_00[it]<<std::setw(9)<<Q_p1[it] 
       <<std::setw(9)<<Q_p2[it]<<" :"<<std::setw(9)<<sum<<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(9)<<qm2<<std::setw(9)<<qm1 
    <<std::setw(9)<<q00<<std::setw(9)<<qp1 
    <<std::setw(9)<<qp2<<" :"<<std::setw(9)<<qt<<std::endl;

//-----------put through amplifiers associated with crosstalks-------------
//  storage arrays for every 0.5 nsec for amplified pulses
    double Pulse[5][100];
 
    double delta_q;
    for (int i=0; i<5; i++) {for (int j=0; j<100; j++) Pulse[i][j]=0.;}

    for (int istrip = 0 ; istrip < 5 ; istrip++ ){
        for (int it = 0 ; it < 50 ; it++ ) {
          switch(istrip){
            case 0: delta_q = Q_m2[it];  break;
            case 1: delta_q = Q_m1[it];  break;
            case 2: delta_q = Q_00[it];  break;
            case 3: delta_q = Q_p1[it];  break;
            case 4: delta_q = Q_p2[it];  break;
          }
          if (delta_q == 0.) continue;
          double time0 = (it + 0.5) * 0.5;
          for (int jt = 1 ; jt < 100 ; jt++ ) {
            double timeAmp = jt * 0.5  ;
            double time = time0 + timeAmp;
            int itime = time / 0.5;
            if( itime >= 100) break;
            double response =  Amp_response(timeAmp)  * delta_q;
            double crosstalk = Amp_crosstalk(timeAmp) * delta_q;
            Pulse[istrip][itime] += response;
            if( istrip > 0) Pulse[istrip-1][itime] += crosstalk;
            if( istrip < 3) Pulse[istrip+1][itime] += crosstalk;
          } 
        }
    }
//   print out Pulse
    std::cout<< std::setprecision(4)<<std::setiosflags(std::ios::fixed)<<
    "=======================amplified pulse================================"
    <<std::endl;
    std::cout<<
    "   time :        Pulse[0] Pulse[1] Pulse[2] Pulse[3] Pulse[4]:Total"
    <<std::endl;
    for (int it = 0; it< 100; it++) {
       double time = (it +1) * 0.5 ;
       double sum =0.;
       for (int i=0;i<4;i++) sum += Pulse[i][it];
       std::cout <<"  "<<time<<" ns : "
       <<std::setw(9)<<Pulse[0][it]<<std::setw(9)<<Pulse[1][it] 
       <<std::setw(9)<<Pulse[2][it]<<std::setw(9)<<Pulse[3][it] 
       <<std::setw(9)<<Pulse[4][it]<<" :"<<std::setw(9)<<sum<<std::endl;
    }
    std::cout<<
    "=======================end of event==================================="
    <<std::endl;

//------------- end of one event loop -------------------------
   return StatusCode::SUCCESS;
}

/////////////////////////////////////////////////////////////////////////

StatusCode MySCTcharge::finalize() {

  std::cout<<"----- Final: number of events = "<<m_eventNumber<<std::endl;
//
return StatusCode::SUCCESS;
}

/////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
//  holeTransport
//---------------------------------------------------------------------
StatusCode MySCTcharge::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 * RandGauss::shoot(m_rndmEngine);
     x += diffusion * RandGauss::shoot(m_rndmEngine);
     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 StatusCode::SUCCESS;
}

//---------------------------------------------------------------------
//  electronTransport
//---------------------------------------------------------------------
StatusCode MySCTcharge::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 * RandGauss::shoot(m_rndmEngine);
     x += diffusion * RandGauss::shoot(m_rndmEngine);
     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 StatusCode::SUCCESS;
}

//----------------------------------------------------------------------
// Electronique response iby Jan Kaplon's form
//     Ic = 220 uA   non-irradiated   signal
//----------------------------------------------------------------------
double MySCTcharge::Amp_response( double t)  {
  if ( t > 0.0) {
  double A=0.31123, B=4.6855, C=0.134685, D=0.0879467, E=8.77451,
         F=0.104167, G=9.39697, H=0.527247, I=0.25;
  double y= 2*(A*cos(D*t)-B*sin(D*t))*exp(-C*t)+E*exp(-F*t)+(-G-H*t)*exp(-I*t);
  y *= 1./0.338;
  return y;
  } else return 0. ;
}

//----------------------------------------------------------------------
// Electronique response iby Jan Kaplon's form
//     Ic = 220 uA   non-irradiated   signal
//----------------------------------------------------------------------
double MySCTcharge::Amp_crosstalk( double t ) {
   if(t > 0.0) {
   double A=0.176508, B=0.538337, C=0.196667, D=2.96078, E=0.134685, 
          F=0.0879467, G=20.1466, H=0.124127, I=8.77451, J=0.104167, 
          K=11.5889, L=0.831773, M=0.25;
   double y = A*exp(-B*t) -2*(C*cos(F*t)-D*sin(F*t))*exp(-E*t)-G*exp(-H*t)
              +I*exp(-J*t)+(K+L*t)*exp(-M*t);
   y *= 1./0.338;
   return y;
   } else return 0.;
}

//---------------------------------------------------------------
//      parameters for electron transport
//---------------------------------------------------------------
bool MySCTcharge::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 MySCTcharge::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 MySCTcharge::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 MySCTcharge::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;
}



//------------------------------------------------------------
//  initialization of e/h transport programme
//-----------------------------------------------------------
void MySCTcharge::initTransport() {

//------------------------ initialize subfunctions ------------
   init_mud_h(m_T) ;
   init_mud_e(m_T) ;
   initExEyArray();
   initPotentialValue();
//----------------------------------------------------------------
   m_kB = 1.38E-23;       // [m^2*kg/s^2/K]
   m_e = 1.602E-19;       // [Coulomb]
//------------ find delepletion deph for model=0 and 1 -------------
   std::cout<<"  model= "<< m_model<<" VB= "<< m_VB <<"  B= "<< m_B << std::endl;
   m_depletion_depth = m_bulk_depth;
   if (m_VB < m_VD) m_depletion_depth = sqrt(m_VB/m_VD) * m_bulk_depth;
// ----------- find for the case of model = 2----------------
   if (m_model == 2) {
     double Ex, Ey;
     for ( double y = 0.; y < 0.0285 ; y += 0.00025 ) {
        EField(0., y, Ex, Ey);   
        if ( Ey > 0.1 ) continue;
        m_depletion_depth = m_bulk_depth - y;
     }
   }
   m_y_origin_min = m_bulk_depth - m_depletion_depth;

   std::cout<<"------ initialization of e-h transport ------"<<std::endl;
   std::cout<<"  m_bulk_depth      = "<< m_bulk_depth <<" [cm]"
            <<std::endl;
   std::cout<<"  m_depletion_depth = "<< m_depletion_depth <<" [cm]"
            <<std::endl;
   std::cout<<"  m_y_origin _min   = "<< m_y_origin_min << " [cm]"
            <<std::endl;

//----for the case of uniform field model (SCT digitization model)------
   double Emean = m_VB / m_depletion_depth;
   m_driftMobility  = mud_h(Emean);
   double r_h = 0.72 - 0.0005*(m_T-273.15);
   m_tanLA = r_h * m_driftMobility * m_B * 1.E-4;
   m_diffusion = m_kB * m_T * m_driftMobility/ m_e;
   std::cout<<"----parameters for old SCT digitization model ----"<< std::endl;
   std::cout<<"  mean E field      = "<< Emean << " [KV/cm] "<< std::endl;
   std::cout<<"  driftMobility     = "<< m_driftMobility<<" [ ]"<< std::endl;
   std::cout<<"  diffusion D       = "<< m_diffusion << " [ ]"<<std::endl;
   std::cout<<"  tan(LA)           = "<< m_tanLA << std::endl;
   std::cout<<"  LA                = "<< atan(m_tanLA)*180./3.141598 
            <<" [degree]"<<std::endl;
   std::cout<<"-----------------------------------------------"<<std::endl;
//--------------------------------------------------------------------

}

//--------------------------------------------------------------
//   drift mobility for electrons
//--------------------------------------------------------------
double MySCTcharge::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 MySCTcharge::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 MySCTcharge::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 MySCTcharge::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 ;
}

//----------------------------------------------------------------------
// Initialize the Random Engine 
//----------------------------------------------------------------------
StatusCode MySCTcharge::initRandomEngine() {
   m_rndmEngine = m_rndmSvc->GetEngine(m_rndmEngineName) ;
   return StatusCode::SUCCESS;
 }

/////////////////////////////////////////////////////////////////////
double ExValue150(int& ix, int& iy), EyValue150(int& ix, int& iy);

void MySCTcharge::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 MySCTcharge::initPotentialValue() {
  for (int ix=0 ; ix <81 ; ix++ ){
     for (int iy=0 ; iy<115 ; iy++ ) {
      m_PotentialValue[ix][iy] = GetPotentialValue(ix,iy);
     }
  }
}