// simulation of induced current model iusing e-h transport // 2010.9.30 T. Kondo(KEK) // y = 0 for HV plane, y = 0.0285 [cm] for strip plane // #include "SCT1/MySCTphi_v2.h" #include "CLHEP/Random/RandFlat.h" #include "CLHEP/Random/RandGauss.h" #include //////////////////////////////////////////////////////// MySCTphi::MySCTphi(const std::string& name, ISvcLocator* pSvcLocator) : AthAlgorithm(name, pSvcLocator), m_rndmSvc("AtRndmGenSvc",name), m_rndmEngineName("MySCTphi") { //---------------setting basic parameters--------------------------- declareProperty("EfieldModel", m_model = 2); // model=2 for FEM field model declareProperty("IsSCTDigiModel", m_isSCTDigiModel = false); declareProperty("TimeOfThreshold", m_timeOfThreshold = 25.0); declareProperty("IsNewPulseShape", m_isNewPulseShape = false); 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("X_origin_min", m_x_origin_min = 0.0000); // [cm] declareProperty("X_origin_max", m_x_origin_max = 0.0080); // [cm] declareProperty("Track_Phi_min", m_phi_min = -16.); declareProperty("Track_Phi_max", m_phi_max = 8.); declareProperty("Track_Eta_min", m_eta_min = -1.); declareProperty("Track_Eta_max", m_eta_max = 1.); declareProperty("Threshold", m_threshold = 1.); // thresold [FC] declareProperty("Ydivision", m_ydivision = 100); // how many loops over y direction declareProperty("eh_pairs", m_eh_pairs = 108.0); // per microd declareProperty("IsLandau", m_isLandau = false); // flag for Landau distr 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] declareProperty("CrossFactor2sides",m_CrossFactor2sides = 0.1); declareProperty("CrossFactorBack",m_CrossFactorBack = 0.07); declareProperty("PeakTime",m_PeakTime = 21); } //////////////////////////////////////////////////////////////// MySCTphi::~MySCTphi() {} ////////////////////////////////////////////////////////////////////////////// StatusCode MySCTphi::initialize() { StatusCode sc=StatusCode::SUCCESS; //////////////////////////////////////////////////////////// std::cout<<"EfieldModel \t"<< m_model << "\t(default = 2)"<regHist("/AANT/hitStrip",m_h_hitStrip); m_h_phi = new TProfile("phi_dist","phi_distrib",48, -16., 8., 0., 20.); m_sc = m_thistSvc->regHist("/AANT/phiDist",m_h_phi); return sc; } //-------------------------------------------------------- // //-------------------------------------------------------- StatusCode MySCTphi::execute() { m_event_number++; //--------------- dE/dX fluctuation ------------------------------ double Landau = 108.; if(m_isLandau) { for (int j=0; j<100; j++) { Landau = gRandom->Landau(72., 11.8) ; if (Landau < 500.) break; } } double eh_pairs = m_eh_pairs * Landau/108. ; //------------------particle track------------------------------ double radian = 180./3.141593; double width = m_phi_max - m_phi_min; double phi = m_phi_min + width * RandFlat::shoot(m_rndmEngine); width = m_eta_max - m_eta_min; double eta = m_eta_min + width * RandFlat::shoot(m_rndmEngine); double z_theta = 2.* atan(exp(-eta)); width = m_x_origin_max - m_x_origin_min; double x_origin = m_x_origin_min + width * RandFlat::shoot(m_rndmEngine); double charge = m_bulk_depth * 1.E4 * eh_pairs *1.6E-19/1.E-15; double efactor = ( charge / cos(phi/radian) / sin(z_theta) ) /m_ydivision ; //-------- first 50 events ---------------------------- if(m_event_number < 50 ) std::cout << "ev#="<=0 && jstrip < 21) qstrip[jstrip][it] += Q_m2[it]; // m1 jstrip++; if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_m1[it]; // 00 jstrip++; if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_00[it]; // p1 jstrip++; if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_p1[it]; // p2 jstrip++; if (jstrip >=0 && jstrip < 21) qstrip[jstrip][it] += Q_p2[it]; } // --- end of time bin loop } // -----------------end of one track loop --------------- //-------------------amplify the charges ---------------------- double pulse[21][3]; // pulse height at m_timeOfThreshold -25,0,+25 ns for (int i=0 ; i<21; i++) { for(int j=0; j<3 ;j++) pulse[i][j] = 0.; } double timeOfThreshold[3]; for (int i=0 ; i<3 ;i++) timeOfThreshold[i] = m_timeOfThreshold + (i-1)*25.0; for (int istrip = 0; istrip < 21 ; istrip++) { for (int it = 0 ; it < 50 ; it++) { double delta_q = qstrip[istrip][it] * efactor; if (delta_q == 0.) continue; double time0 = (it + 0.5) * 0.5; for (int j=0; j<3 ;j++) { double timeAmp = timeOfThreshold[j] - time0; if ( timeAmp <= 0.) continue; double response = Amp_response(timeAmp) * delta_q; double crosstalk = Amp_crosstalk(timeAmp) * delta_q; pulse[istrip][j] += response; if( istrip > 0) pulse[istrip-1][j] += crosstalk; if( istrip < 19) pulse[istrip+1][j] += crosstalk; } } } if(m_coutLevel>1) { std::cout<<"tbin\ttime\tpulse[8]\tpulse[9]\tpulse[10]\tpulse[11]" <Fill(istrip-10+0.5); } } // ----- find maximum cluster size --------------- double max_clsize = 0.; bool previous = false; double cluster_size = 0.; for(int istrip=0 ;istrip<21;istrip++ ) { if( tbin[istrip][1] == 1 ) { if (!previous) { cluster_size = 1.; previous = true; } else cluster_size += 1.; } else { previous = false; cluster_size = 0.; } max_clsize = (cluster_size>max_clsize) ? cluster_size : max_clsize; } if( m_coutLevel>0 || fmod(m_event_number,2000) == 0) std::cout<GetEngine(m_rndmEngineName) ; // if (m_rndmEngine==0) { // ATH_MSG_ERROR ( "Could not find RndmEngine : " << m_rndmEngineName ) ; // return StatusCode::FAILURE ; // } // ATH_MSG_INFO ( "Get random number engine : <" << m_rndmEngineName << ">" ) ; return StatusCode::SUCCESS; } //---------------------------------------------------------------------- // Electronique response // new Amp : by Jan Kaplon's form with Ic = 220 uA (non-irradiated ) // old Amp : CR-RC^3 of the charge diode //---------------------------------------------------------------------- double MySCTphi::Amp_response( double time) { double y = 0.; if ( time > 0.0) { if ( m_isNewPulseShape ) { double t = time; 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; 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; } else { double tC = time / (m_PeakTime/3.0) ; y = tC*tC*tC*exp(-tC) * m_NormConstCentral ; } } return y ; } //---------------------------------------------------------------------- // differenciated and scaled pulse on the neighbour strip! //---------------------------------------------------------------------- double MySCTphi::Amp_crosstalk( double time ) { double y = 0.; if(time > 0.0) { if ( m_isNewPulseShape ) { double t = time; 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; 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; } else { double tC = time / (m_PeakTime/3.0) ; y = tC*tC*exp(-tC)*(3.0-tC) * m_NormConstNeigh ; } } return y; } //---------------------------------------------------------------------- // Initialize Amplifier //---------------------------------------------------------------------- void MySCTphi::init_Amp() { m_NormConstCentral = (exp(3.0)/27.0)*(1.0-m_CrossFactor2sides)*(1.0-m_CrossFactorBack); 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 ; } ///////////////////////////////////////////////////////////////////////// //--------------------------------------------------------------------- // holeTransport //--------------------------------------------------------------------- StatusCode MySCTphi::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="<(" <(" <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 = kB * m_T * mu_e/ e; return true; } else return false; } //--------------------------------------------------------------- // parameters for hole transport //--------------------------------------------------------------- bool MySCTphi::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 MySCTphi::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_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 MySCTphi::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 MySCTphi::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 ; } ///////////////////////////////////////////////////////////////////// //------------------------------------------------------------ // initialization of e/h transport programme //----------------------------------------------------------- void MySCTphi::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 ------"<