// 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 "SCT2/MySCTtscan_v1.h" #include "CLHEP/Random/RandFlat.h" #include "CLHEP/Random/RandGauss.h" #include //////////////////////////////////////////////////////// MySCTtscan::MySCTtscan(const std::string& name, ISvcLocator* pSvcLocator) : AthAlgorithm(name, pSvcLocator), m_rndmSvc("AtRndmGenSvc",name), m_rndmEngineName("MySCTtscan") { //---------------setting basic parameters--------------------------- declareProperty("EfieldModel", m_model = 2); // model=2 for FEM field model declareProperty("IsSCTDigiModel", m_isSCTDigiModel = false); declareProperty("TimeOffset", m_timeOffset = 25.0); declareProperty("Noise", m_noise = 0.); // [number of electrons]^ declareProperty("TimeOfThreshold_min", m_timeOfThreshold_min = -25.0); declareProperty("TimeOfThreshold_max", m_timeOfThreshold_max = 75.0); declareProperty("ChargeTimeStep", m_chargeTimeStep = 0.5 ); // declareProperty("IsNewPulseShape", m_isNewPulseShape = true); declareProperty("IsLevelSensing", m_isLevelSensing = true); //false is edge-sensing declareProperty("Felectron", m_felectron =0.00) ; declareProperty("Fpion", m_fpion = 1.00) ; declareProperty("Fkaon", m_fkaon =0.00) ; declareProperty("Fproton", m_fproton =0.00) ; declareProperty("StripNo_min", m_strip_min = 0) ; declareProperty("StripNo_max", m_strip_max = 768) ; declareProperty("LocPhi_min" , m_locf_min = -30. ); declareProperty("LocPhi_max" , m_locf_max = 20. ); declareProperty("Polarity" , m_polarity = 0 ); // 0 for + and - declareProperty("Pt_min" , m_pt_min = 0.5 ) ; declareProperty("Pt_max" , m_pt_max = 100. ) ; declareProperty("PtPower" , m_ptPower = -7.0) ; declareProperty("PtPowerForSampling" , m_ptPowerForSampling = -2.0) ; declareProperty("BarrelLayer" , m_layer = 0 ) ; declareProperty("Zmodule_min" , m_zmodule_min = 6.0288-6.2 ) ; declareProperty("Zmodule_max" , m_zmodule_max = 6.0288+6.2 ) ; 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("StripTilt", m_stripTilt = 1);//0=no tilt,1=50-50,2=all tilt declareProperty("Threshold", m_threshold = 1.); // thresold [FC] declareProperty("CrossTalk", m_crossTalk = 1.); // cross talk factor declareProperty("Ydivision", m_ydivision = 100); // how many loops over y direction declareProperty("eh_pairs", m_eh_pairs = 108.0); // per microd declareProperty("IsdEdX", m_isdEdX = false); // flag for Landau distr 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); } //////////////////////////////////////////////////////////////// MySCTtscan::~MySCTtscan() {} ////////////////////////////////////////////////////////////////////////////// StatusCode MySCTtscan::initialize() { StatusCode sc=StatusCode::SUCCESS; m_event_number = 0; m_sc=StatusCode::SUCCESS; ATH_MSG_INFO ("initialize()"); //+++ Get the random generator engine m_sc = initRandomEngine(); if(m_sc.isFailure()) return StatusCode::FAILURE ; // initLocPhi() ; // initTransport(); init_Amp(); m_noise_FC = m_noise * 1.602E-19/1.E-15; /** get a handle on the NTuple and histogramming service */ m_sc = service("THistSvc", m_thistSvc); if (m_sc.isFailure()) { ATH_MSG_ERROR ( "Unable to retrieve pointer to THistSvc"); return m_sc; } //////////////////////////////////////////////////////////// std::cout<<"EfieldModel \t"<< m_model << "\t(default = 2)"<regHist("/AANT/hitStrip",m_h_hitStrip); m_h_pt = new TH1F("pt", "pt",100, -2., 2.); m_sc = m_thistSvc->regHist("/AANT/pt",m_h_pt); m_h_locf = new TH1F("locf", "locf",90, -90., 90.); m_sc = m_thistSvc->regHist("/AANT/locf",m_h_locf); double t0 = m_timeOfThreshold_min - 2.5 - m_timeOffset; double t1 = m_timeOfThreshold_max - 2.5 - m_timeOffset; m_h_tscan_0 = new TProfile("T1XY","tbin_0",20, t0, t1, 0., 20.); m_sc = m_thistSvc->regHist("/AANT/T1XY",m_h_tscan_0); m_h_tscan_1 = new TProfile("T01X","tbin_1",20, t0, t1, 0., 20.); m_sc = m_thistSvc->regHist("/AANT/T01X",m_h_tscan_1); m_h_tscan_2 = new TProfile("T001","tbin_2",20, t0, t1, 0., 20.); m_sc = m_thistSvc->regHist("/AANT/T001",m_h_tscan_2); m_h_tscan_all = new TProfile("TALL","tbin_all",20, t0,t1, 0., 20.); m_sc = m_thistSvc->regHist("/AANT/TALL",m_h_tscan_all); //-----get pt distribution from ifile pt.root ------------------- m_h_ptin = new TH1F("ptin","ptin",100, 0., 2.15); m_sc = m_thistSvc->regHist("/AANT/ptin",m_h_ptin); ifstream fin("pt.dat"); int idummy; double x,y; while (fin>>idummy>>x>>y ) { std::cout << idummy <<" "<< x <<" "<< y << std::endl; m_h_ptin->Fill(x/1000.,y); } fin.close(); //--------------------------------------------------- m_ptArea = (pow(m_pt_max,1.+m_ptPowerForSampling) - pow(m_pt_min,m_ptPowerForSampling+1.)) /(m_ptPowerForSampling+1.); return sc; } //-------------------------------------------------------- // //-------------------------------------------------------- StatusCode MySCTtscan::execute() { m_event_number++; //-----loop over m_timeOfThreshold_min to m_timeOfThreshold_max ----------- double timeOfThreshold[3]; for (int itot=0 ; itot < 100 ; itot++) { m_timeOfThreshold = m_timeOfThreshold_min + itot * 5.0 ; double trueTime = m_timeOfThreshold - m_timeOffset; if( m_timeOfThreshold >= m_timeOfThreshold_max ) break; if( m_coutLevel > 1) std::cout<<"m_timeOfThreshold = "<GetRandom(); if(m_polarity<0 || ((m_polarity==0)&&(RandFlat::shoot(m_rndmEngine)< 0.5))) pt = -pt; m_h_pt->Fill(pt); //----------------select z on module-------------------- double z_module = m_zmodule_min + (m_zmodule_max-m_zmodule_min) * RandFlat::shoot(m_rndmEngine); double z_IP = 7.0 * RandGauss::shoot(m_rndmEngine); double z_theta = atan2(m_radius, z_module - z_IP); double eta = -log(tan(z_theta/2.)); // selection of stereo angle bool isStripTilt = false; if(m_stripTilt > 0) { if(m_stripTilt ==2 ) isStripTilt = true; else if( RandFlat::shoot(m_rndmEngine) > 0.5) isStripTilt = true; } double stereo = 0.; if (isStripTilt) stereo = 0.04; //----- select emmision angle in phi in the xy plane--> x_origin, locf double locf=0.; double x_origin = 0.; for (int iphi0 = 0; iphi0 < 5000 ; iphi0++) { if (iphi0==4999) { std::cout<<"***** CANNOT FIND GOOD PHI0 IN 499999 TIMES****"<0) phi0 = (m_phi0_max_pos - m_phi0_min_pos)*ran + m_phi0_min_pos; if(pt<0) phi0 = (m_phi0_max_neg - m_phi0_min_neg)*ran + m_phi0_min_neg; phi0 = phi0/radian; locf = locPhi(stereo, pt, phi0, z_theta, x_origin) * radian; if (locf > 900. ) continue; // check if locf is in the range specified if( m_coutLevel>0) std::cout <<"n,stereo,locf,x_origin= "< m_locf_max) continue; // limit events within the module --- if (x_origin > m_x_origin_min && x_origin < m_x_origin_max) break; } m_h_locf->Fill(locf); //-------------calculate TOF ---------------------------- //----------- get p from pt ----------------- double z_factor = 1./ sin(z_theta); double p = pt * z_factor; // int id_particle = 1; // 0(e), 1(pion), 2(K), 3(p) double fp = RandFlat::shoot(m_rndmEngine); if ( fp > m_fpion) id_particle = 0; if ( fp > m_fpion + m_felectron) id_particle = 2; if ( fp > m_fpion + m_felectron + m_fkaon ) id_particle = 3; double dedx_factor = 1.; if(m_isdEdX) dedx_factor = dEdX(id_particle, p); //--------------------------------------------- double qstrip[21][50]; for (int i=0 ; i<21; i++) { for(int j=0; j<50 ;j++) qstrip[i][j] = 0.;} //--------------- dE/dX fluctuation ------------------------------ double LandauFactor = 1.; if(m_isLandau) { double Landau; for (int j=0; j<100; j++) { Landau = gRandom->Landau(72., 11.8) ; if (Landau < 500.) break; } LandauFactor = Landau / 108.; } double eh_pairs = m_eh_pairs * LandauFactor * dedx_factor ; double charge = m_bulk_depth * 1.E4 * eh_pairs *1.6E-19/1.E-15; double efactor = ( charge * z_factor / cos(locf/radian) ) /m_ydivision ; if(m_coutLevel>0) std::cout<<"event,id,pt,eta,locf,z_factor,dedx,x_origin,eh_pairs=" < "<1) { std::cout<<"tbin\ttime\tpulse[8]\tpulse[9]\tpulse[10]\tpulse[11]" <Fill(istrip-10+0.5); } if( m_coutLevel>0 || fmod(m_event_number,2000) == 0) std::cout<0 && tbin[istrip][itbin-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( max_clsize > 0) { if( itbin ==0 ) m_h_tscan_0 -> Fill(trueTime,max_clsize,1); if( itbin ==1 ) m_h_tscan_1 -> Fill(trueTime,max_clsize,1); if( itbin ==2 ) m_h_tscan_2 -> Fill(trueTime,max_clsize,1); } } // ---- find cluster size for 1XX+X1X+XX1 --------------- double max_clsize = 0.; bool previous = false; double cluster_size = 0.; for(int istrip=0 ;istrip<21;istrip++ ) { if(tbin[istrip][0]+tbin[istrip][1]+tbin[istrip][2] > 0) { 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( max_clsize > 0) m_h_tscan_all -> Fill(trueTime,max_clsize,1); // } // end of loop for m_timeOverThreshold scan return StatusCode::SUCCESS; } ///////////////////////////////////////////////////////////////////////// StatusCode MySCTtscan::finalize() { MsgStream mLog(msgSvc(), name()); mLog << MSG::INFO << "finalize()" << endreq; std::cout <<"Final Event Number = "<GetXaxis()->GetBinCenter(i+1); double n_0 = m_h_tscan_0->GetBinEntries(i+1); double n_1 = m_h_tscan_1->GetBinEntries(i+1); double n_2 = m_h_tscan_2->GetBinEntries(i+1); double n_all = m_h_tscan_all->GetBinEntries(i+1); double mean_0 = m_h_tscan_0->GetBinContent(i+1); double mean_1 = m_h_tscan_1->GetBinContent(i+1); double mean_2 = m_h_tscan_2->GetBinContent(i+1); double mean_all = m_h_tscan_all->GetBinContent(i+1); double smean_0 = m_h_tscan_0->GetBinError(i+1); double smean_1 = m_h_tscan_1->GetBinError(i+1); double smean_2 = m_h_tscan_2->GetBinError(i+1); double smean_all = m_h_tscan_all->GetBinError(i+1); std::cout<<"tot= "<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 MySCTtscan::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 MySCTtscan::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 MySCTtscan::init_Amp() { //m_NormConstCentral = (exp(3.0)/27.0)*(1.0-m_CrossFactor2sides)*(1.0-m_CrossFactorBack); m_NormConstCentral = (exp(3.0)/27.0); 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); m_NormConstNeigh *= (m_CrossFactor2sides/2.0); return ; } ///////////////////////////////////////////////////////////////////////// //--------------------------------------------------------------------- // holeTransport //--------------------------------------------------------------------- StatusCode MySCTtscan::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 MySCTtscan::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 MySCTtscan::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 MySCTtscan::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 MySCTtscan::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 MySCTtscan::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 ------"<