//----------------------------------------------- // 2016.4.3 for ICM (induced charge model) by Taka Kondo (KEK) //----------------------------------------------- #include "GetExEy_150V.cxx" #include "GetPotential.cxx" void initialize(int set_m_EfieldModel, double set_m_VD, double set_m_VB) ; 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() ; double Amp_response( double t) ; double Amp_crosstalk( double t ) ; void setupCanvas(const int NX, const int NY) ; TCanvas* canvas; TPad* canv[5]; //------parameters given exterbnally by jobOptions ------------------ int m_EfieldModel ; 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 = 1.38E-23; // [m^2*kg/s^2/K] double m_e = 1.602E-19; // [Coulomb] double m_vs_e; double m_Ec_e ; double m_vs_h; double m_Ec_h ; double m_driftMobility; // [cm**2/V/s] double m_diffusion; // [cm**2/s] double m_beta_e ; double m_beta_h ; double m_theta; double m_tanLA; //---------parameters for amplifier ---------------------------- double m_PeakTime = 21.; // [ns] double m_CrossFactor2sides = 0.10 ; // 10% cross talk (5% each) double m_NormConstCentral = (exp(3.0)/27.0); // SCT Digitization model double m_NormConstNeigh = exp(3.0-sqrt(3.0))/(6*(2.0*sqrt(3.0)-3.0)) * (m_CrossFactor2sides/2.0) ; // normalize //---------- arrays of FEM analysis ----------------------------------- double m_PotentialValue[81][115]; double m_ExValue150[17][115]; double m_EyValue150[17][115]; //-----------------char string for sprintf------------------------------- char m_cid[200]; //-------histogramme ---------------------------------- void initialize(int set_m_EfieldModel, double set_m_VD, double set_m_VB) { //---------------setting basic parameters--------------------------- m_isSCTDigiModel = false ; // true for SCT diditization model (as of 2010)a // false for ICM (induced charge model) m_EfieldModel = set_m_EfieldModel ; // 0 (uniform E-field), 1 (flat diode model), 2 (FEM solusions) m_VD = set_m_VD ; // full depletion voltage [Volt] m_VB = set_m_VB ; // applied bias voltage [Volt] ////------------------usually fixed------------------------------------- m_bulk_depth = 0.0285 ; // in [cm] m_strip_pitch = 0.0080 ; // in [cm] m_T = 273.15; m_B = -2.0 ; // [Tesla] //---------------------------------------------------------------- m_transportTimeStep = 0.25; // one step side in time [nsec] m_transportTimeMax = 25.0; // maximun tracing time [nsec] m_coutLevel = 0 ; // specify printout level //------------ find delepletion deph for model=0 and 1 ------------- m_depletion_depth = m_bulk_depth; if (m_VB < m_VD) m_depletion_depth = sqrt(m_VB/m_VD) * m_bulk_depth; //-------------------------------------------------------- std::cout<<"----------------------------------------------"<(" <(" <xfinal,istrip,dx=" <"<t_current="< 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 = m_kB * m_T * mu_e/ m_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 = m_kB * m_T * mu_h/ m_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="< 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_EfieldModel ==0 ) { if ( m_bulk_depth - y < m_depletion_depth ) { Ey = m_VB / m_depletion_depth ; } else { Ey = 0.; } return; } //---------- case for flat diode model ------------------------------ if(m_EfieldModel==1) { if(m_VB > m_VD) { //Ey = (m_VB+m_VD)/m_depletion_depth Ey = (m_VB+m_VD)/m_bulk_depth - 2.*m_VD*(m_bulk_depth-y)/(m_bulk_depth*m_bulk_depth); } else { if ( m_bulk_depth - y < m_depletion_depth ) { double Emax = 2.* m_depletion_depth * m_VD / (m_bulk_depth*m_bulk_depth); Ey = Emax*(1-(m_bulk_depth-y)/m_depletion_depth); } else { Ey = 0.; } } 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 ; } //---------------------------------------------------------------------- //---------------------------------------------------------------------- double Amp_response( double t ) { double y = 0.; if ( t > 0.0) { double tC = t / (m_PeakTime/3.0) ; y = tC*tC*tC*exp(-tC) * m_NormConstCentral ; //------------------------------------------------- // Electronique response by Jan Kaplon's form // Ic = 220 uA non-irradiated signal // 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; } //---------------------------------------------------------------------- double Amp_crosstalk( double t ) { double y = 0.; if(t > 0.0) { double tC = t / (m_PeakTime/3.0) ; y = tC*tC*exp(-tC)*(3.0-tC) * m_NormConstNeigh ; } //------------------------------------------------- // Electronique response by Jan Kaplon's form // Ic = 220 uA non-irradiated signal // 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; }