//----------------------------------------------- // 2017.7.24 update for the case of negative VD (type-P) // 2017.7.17 updated // 2016.4.3 for ICM (induced charge model) by Taka Kondo (KEK) //----------------------------------------------- #include "GetExEy_150V.cxx" #include "GetPotential.cxx" 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() ; void init_Amp(); double Amp_response( double t) ; double Amp_crosstalk( double t ) ; //------parameters given externally by jobOptions ------------------ int m_EfieldModel = 1 ; // 0(uniform E), 1(flat diode model), 2 (FEM solusions) bool m_isSCTDigiModel = false ; // true (SCTDigi model) false (Ramo model) bool m_isNewCrossTalk = true ; // true(newCrossTalk) false(SCTDigi crosstalk) double m_VD = 70. ; // full depletion voltage [Volt] negative for type-P double m_VB = 150. ; // applied bias voltage [Volt] double m_T = 273.15; double m_B = -2.0 ; // [Tesla] double m_transportTimeStep = 0.50; // one step side in time [nsec] double m_transportTimeMax = 50.0; // maximun tracing time [nsec] double m_x0; double m_y0; int m_coutLevel = 0 ; int m_eventNumber = 0 ; //------parameters mostly fixed but can be changed externally ------------ double m_bulk_depth = 0.0285 ; // in [cm] double m_strip_pitch = 0.0080 ; // in [cm] 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 particle trajectories ----------------------- double m_x_origin_min = 0.0000 ; double m_x_origin_max = 0.0080 ; double m_phi_min = -16.; double m_phi_max = 8. ; double m_eta_min = -1. ; double m_eta_max = 1. ; double m_threshold = 1. ; // threshold [fC] double m_eh_pairs = 73.0 ; // most probable number of eh_pairs per micron int m_ydivision = 100; // how many loops over y direction bool m_isLandau = false ; //---------parameters for amplifier ---------------------------- double m_PeakTime = 21.; // [ns] double m_NormConstCentral ; double m_NormConstNeigh ; double m_CrossFactor2sides = 0.10 ; // 10% cross talk (5% each) double m_CrossFactorBack = 0.07; double m_timeOfThreshold = 25. ; //---------parameters for 25ns gate flag and its gate timing -------------- bool m_is25nsGate = false; double m_GateStartTiming = 14.0; // for Ramo model (20ns for SCTDigi) double m_GateWidth = 25.0 ; // 25ns gate width [ns] //---------- 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() { //---------------setting basic parameters--------------------------- //m_isSCTDigiModel = false ; // true for SCT diditization model (as of 2010)a // false for ICM (induced charge model) //------------------usually fixed------------------------------------- if( m_isSCTDigiModel ) { m_timeOfThreshold = 30.0; // 30ns for original model m_GateStartTiming = 20.0 ; // 20ns for SCTDigi model } else { m_timeOfThreshold = 25.0; // 30ns for original model m_GateStartTiming = 14.0 ; // 20ns for SCTDigi model } //------------------------ initialize subfunctions ------------ init_mud_h(m_T) ; init_mud_e(m_T) ; init_Amp(); initExEyArray(); initPotentialValue(); //------------ find delepletion deph for model=0 and 1 ------------- // for type N (before type inversion) if (m_VD >= 0) { m_depletion_depth = m_bulk_depth; if (m_VB < m_VD) m_depletion_depth = sqrt(m_VB/m_VD) * m_bulk_depth; } else { // for type P m_depletion_depth = m_bulk_depth; if ( m_VB <= abs( m_VD) ) m_depletion_depth = 0.; } //-------------------------------------------------------- if (m_coutLevel >=0) { std::cout<<"----------------------------------------------"<=0) { 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<<" [cm**2/V/s]"<< std::endl; std::cout<<" diffusion D = "<< m_diffusion << " [cm**2/s]"<= 0) { std::cout<<"---- parameters for electron transport -----"<= 0) { std::cout<<"----parameters for hole transport -----"<1) std::cout<<"holeTransport : x,y=("<< x0*1.e4<<","<(" <1) std::cout<<"elecTransport : x,y=("<< x0*1.e4<<","<(" < 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 of induced charge using Weighting (Ramo) function //------------------------------------------------------------------- double induced (int istrip, double x, double y) { // x and y are the location of charge (e or hole) // induced charge 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 && m_depletion_depth >0. ) { Ey = m_VB / m_depletion_depth ; } else { Ey = 0.; } return; } //---------- case for flat diode model ------------------------------ if(m_EfieldModel==1) { if(m_VB > abs(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 && m_depletion_depth >0.) { //double Emax = 2.* m_depletion_depth * m_VD / double Emax = 2.* m_depletion_depth * abs(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 (for SCTDigi model) //---------------------------------------------------------------------- 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 (for SCTDigi model) //---------------------------------------------------------------------- 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 ; } //---------------------------------------------------------------------- // Initialize Amplifier //---------------------------------------------------------------------- void init_Amp() { m_NormConstCentral = (exp(3.0)/27.0) ; m_NormConstCentral *= (1.0-m_CrossFactor2sides)*(1.0-m_CrossFactorBack); if( m_isNewCrossTalk ) m_NormConstNeigh = (exp(3.0)/27.0); else m_NormConstNeigh = exp(3.0-sqrt(3.0))/(6*(2.0*sqrt(3.0)-3.0)); m_NormConstNeigh *= (m_CrossFactor2sides/2.0)*(1.0-m_CrossFactorBack); return ; } //---------------------------------------------------------------------- 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) ; if ( m_isNewCrossTalk ) y = tC*tC*tC*exp(-tC) * m_NormConstNeigh ; else 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; }