diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt b/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
index 5aecd81ce267de181e1ed19ccd4f2d4caf2eb1b8..af369ec9c0118b0e5960eefb3cd46b059d9008eb 100644
--- a/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
@@ -27,6 +27,7 @@ atlas_depends_on_subdirs( PUBLIC
InnerDetector/InDetDetDescr/InDetReadoutGeometry
InnerDetector/InDetDetDescr/SCT_ModuleDistortions
InnerDetector/InDetRawEvent/InDetSimData
+ MagneticField/MagFieldInterfaces
Simulation/Tools/AtlasCLHEP_RandomGenerators )
# External dependencies:
@@ -39,7 +40,7 @@ atlas_add_component( SCT_Digitization
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesSvcLib InDetIdentifier InDetReadoutGeometry InDetSimData AtlasCLHEP_RandomGenerators )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesSvcLib InDetIdentifier InDetReadoutGeometry InDetSimData AtlasCLHEP_RandomGenerators MagFieldInterfaces PathResolver )
# Install files from the package:
atlas_install_headers( SCT_Digitization )
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c481fa99a25bf529a6be08e8c789affb6625098f
--- /dev/null
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx
@@ -0,0 +1,453 @@
+//-----------------------------------------------
+// 2020.8.12 Implementation in Athena by Manabu Togawa
+// 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 "SCT_InducedChargedModel.h"
+
+void SCT_InducedChargedModel::Init(float vdepl, float vbias, IdentifierHash m_hashId, ServiceHandle<ISiliconConditionsSvc> m_siConditionsSvc) {
+
+ ATH_MSG_INFO("--- Induced Charged Model Paramter (Begin) --------");
+
+//---------------setting basic parameters---------------------------
+
+ m_VD = vdepl; // full depletion voltage [Volt] negative for type-P
+ m_VB = vbias; // applied bias voltage [Volt]
+ m_T = m_siConditionsSvc->temperature(m_hashId) + Gaudi::Units::STP_Temperature;
+
+//------------------------ initialize subfunctions ------------
+
+ loadICMParameters();
+
+//------------ 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.;
+ }
+
+//--------------------------------------------------------
+
+ ATH_MSG_INFO(" EfieldModel 0(uniform E), 1(Flat Diode Model), 2 (FEM solusions)\t\t"<< m_EfieldModel );
+ ATH_MSG_INFO(" DepletionVoltage_VD\t"<< m_VD <<" V");
+ ATH_MSG_INFO(" BiasVoltage_VB \t"<< m_VB <<" V");
+ ATH_MSG_INFO(" MagneticField_B \t"<< m_B <<" Tesla");
+ ATH_MSG_INFO(" Temperature_T \t"<< m_T <<" K");
+ ATH_MSG_INFO(" TransportTimeStep \t"<< m_transportTimeStep <<" ns");
+ ATH_MSG_INFO(" TransportTimeMax\t"<< m_transportTimeMax <<" ns");
+ ATH_MSG_INFO(" BulkDepth\t\t"<< m_bulk_depth << " cm");
+ ATH_MSG_INFO(" DepletionDepth\t\t"<< m_depletion_depth << " cm");
+ ATH_MSG_INFO(" StripPitch\t\t"<< m_strip_pitch << " cm");
+ ATH_MSG_INFO("--- Induced Charged Model Paramter (End) ---------");
+
+ return;
+}
+
+//---------------------------------------------------------------------
+// holeTransport
+//---------------------------------------------------------------------
+void SCT_InducedChargedModel::holeTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) {
+
+// transport holes in the bulk
+// T. Kondo, 2010.9.9, 2017.7.20
+// External parameters to be specified
+// m_transportTimeMax [nsec]
+// m_transportTimeStep [nsec]
+// m_bulk_depth [cm]
+// Induced currents are added to every m_transportTimeStep (defailt is 0.5 ns):
+// Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100]
+// means 50 ns storages of charges in each strips.
+
+ double x = x0; // original hole position [cm]
+ double y = y0; // original hole position [cm]
+ bool isInBulk = true;
+ double t_current = 0.;
+ double qstrip[5]; // induced charges on strips due to a hole
+ double vx, vy, D;
+ for (int istrip = -2 ; istrip <= 2 ; istrip++) qstrip[istrip+2] =
+ induced(istrip, x, y); // original induced charge bu hole +e
+ while ( t_current < m_transportTimeMax ) {
+ if ( !isInBulk ) break;
+ if ( !hole( x, y, vx, vy, D, m_hashId, m_siPropertiesSvc)) break ;
+ double delta_y = vy * m_transportTimeStep *1.E-9;
+ y += delta_y;
+ double dt = m_transportTimeStep;
+ if ( y > m_bulk_depth ) {
+ isInBulk = false ; // outside bulk region
+ 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 * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
+ x += diffusion * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
+ if ( y > m_bulk_depth ) {
+ y = m_bulk_depth;
+ isInBulk = false; // outside bulk region
+ }
+
+// get induced current by subtracting induced charges
+ for (int istrip = -2 ; istrip <= 2 ; istrip++) {
+ double qnew = induced( istrip, x, y);
+ int jj = istrip + 2;
+ double dq = qnew - qstrip[jj];
+ qstrip[jj] = qnew ;
+
+ int jt = int( (t_current+0.001) / m_transportTimeStep) ;
+ if (jt < 100) {
+ 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
+
+ return ;
+}
+
+//---------------------------------------------------------------------
+// electronTransport
+//---------------------------------------------------------------------
+void SCT_InducedChargedModel::electronTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc ) {
+
+// transport electrons in the bulk
+// T. Kondo, 2010.9.9, 2017.7.20
+// External parameters to be specified
+// m_transportTimeMax [nsec]
+// m_transportTimeStep [nsec]
+// m_bulk_depth [cm]
+// Induced currents are added to every m_transportTimeStep (defailt is 0.5 ns):
+// Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100]
+// means 50 ns storages of charges in each strips.
+
+ 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 <= 2 ; istrip++)
+ qstrip[istrip+2] = -induced(istrip, x, y);
+ while (t_current < m_transportTimeMax) {
+ if ( !isInBulk ) break;
+ if ( !electron( x, y, vx, vy, D, m_hashId, m_siPropertiesSvc)) 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 * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
+ x += diffusion * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
+ 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 <= 2 ; istrip++) {
+ double qnew = -induced( istrip, x, y);
+ int jj = istrip + 2;
+ double dq = qnew - qstrip[jj];
+ qstrip[jj] = qnew ;
+
+ int jt = int( (t_current + 0.001) / m_transportTimeStep);
+ if (jt< 100) {
+ 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 electron tracing
+
+ return ;
+}
+
+//---------------------------------------------------------------
+// parameters for electron transport
+//---------------------------------------------------------------
+bool SCT_InducedChargedModel::electron(double x_e, double y_e, double &vx_e, double &vy_e, double &D_e, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) {
+
+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 = m_siPropertiesSvc->getSiProperties(m_hashId).calcElectronDriftMobility(m_T,E*CLHEP::volt/CLHEP::cm);
+ mu_e *= (CLHEP::volt/CLHEP::cm)/(CLHEP::cm/CLHEP::s);
+ v_e = mu_e * E;
+ r_e = 0;
+ r_e = 1.13+0.0008*(m_T-Gaudi::Units::STP_Temperature);
+ 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 SCT_InducedChargedModel::hole(double x_h, double y_h, double &vx_h, double &vy_h, double &D_h, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) {
+
+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 = m_siPropertiesSvc->getSiProperties(m_hashId).calcHoleDriftMobility(m_T,E*CLHEP::volt/CLHEP::cm);
+ mu_h *= (CLHEP::volt/CLHEP::cm)/(CLHEP::cm/CLHEP::s);
+ v_h = mu_h * E;
+ r_h = 0;
+ r_h = 0.72 - 0.0005*(m_T-Gaudi::Units::STP_Temperature);
+ 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 SCT_InducedChargedModel::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 ;
+
+ return P;
+}
+
+//--------------------------------------------------------------
+// Electric Field Ex and Ey
+//--------------------------------------------------------------
+void SCT_InducedChargedModel::EField( double x, double y, double &Ex, double &Ey ) {
+// m_EfieldModel == 0; uniform E-field model
+// m_EfieldModel == 1; flat diode model
+// m_EfieldModel == 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_EfieldModel==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;
+
+ 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-----------------------
+
+ Ex00 = m_ExValue[ix][iy]; Ex10 = m_ExValue[ix1][iy];
+ Ex01 = m_ExValue[ix][iy1]; Ex11 = m_ExValue[ix1][iy1];
+ Ey00 = m_EyValue[ix][iy]; Ey10 = m_EyValue[ix1][iy];
+ Ey01 = m_EyValue[ix][iy1]; Ey11 = m_EyValue[ix1][iy1];
+
+//---------------- end of data bank search---
+ Ex = Ex00*(1.-fx)*(1.-fy) + Ex10*fx*(1.-fy)
+ + Ex01*(1.-fx)*fy + Ex11*fx*fy ;
+ 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_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_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 * abs(m_VD) /
+ (m_bulk_depth*m_bulk_depth);
+ Ey = Emax*(1-(m_bulk_depth-y)/m_depletion_depth);
+ } else {
+ Ey = 0.;
+ }
+ }
+ return;
+ }
+return;
+}
+
+/////////////////////////////////////////////////////////////////////
+
+void SCT_InducedChargedModel::loadICMParameters(){
+// Loading Ramo potential and Electric field.
+// In strip pitch directions :
+// Ramo potential : 80 divisions (81 points) with 5 um intervals from 40-440 um.
+// Electric field : 16 divisions (17 points) with 2.5 um intervals from 0-40 um.
+// In sensor depth directions (for both potential and Electric field):
+// 114 divisions (115 points) with 2.5 nm intervals for 285 um.
+
+ TFile *hfile = new TFile( PathResolverFindCalibFile("SCT_Digitization/SCT_InducedChargedModel.root").c_str() );
+
+ h_Potential_FDM = static_cast<TH2F *>(hfile->Get("Potential_FDM"));
+ h_Potential_FEM = static_cast<TH2F *>(hfile->Get("Potential_FEM"));
+
+ constexpr float MinBiasV_ICM_FEM = -180.0;
+ constexpr float MaxBiasV_ICM_FEM = 70.0;
+ if ( m_VD<MinBiasV_ICM_FEM || m_VD>MaxBiasV_ICM_FEM ){
+ m_EfieldModel = 1; // Change to FDM
+ ATH_MSG_INFO("Changed to Flat Diode Model since deplettion volage is out of range. (-180 < m_VD < 70 is allow.)");
+ }
+
+ // For Ramo Potential
+ for (int ix=0 ; ix <81 ; ix++ ){
+ for (int iy=0 ; iy<115 ; iy++ ) {
+ if (m_EfieldModel==2){ // FEM
+ m_PotentialValue[ix][iy] = h_Potential_FEM -> GetBinContent(ix+1,iy+1);
+ } else { // FDM
+ m_PotentialValue[ix][iy] = h_Potential_FDM -> GetBinContent(ix+1,iy+1);
+ }
+ }
+ }
+
+ h_Potential_FDM -> Delete();
+ h_Potential_FEM -> Delete();
+
+ if ( m_EfieldModel == 2 ){
+
+ float VFD0[9]= { -180, -150, -120, -90, -60, -30, 0, 30, 70};
+
+ std::vector<std::string> hx_list;
+ std::vector<std::string> hy_list;
+
+ hx_list.push_back("Ex_FEM_M180_0");
+ hx_list.push_back("Ex_FEM_M150_0");
+ hx_list.push_back("Ex_FEM_M120_0");
+ hx_list.push_back("Ex_FEM_M90_0");
+ hx_list.push_back("Ex_FEM_M60_0");
+ hx_list.push_back("Ex_FEM_M30_0");
+ hx_list.push_back("Ex_FEM_0_0");
+ hx_list.push_back("Ex_FEM_30_0");
+ hx_list.push_back("Ex_FEM_70_0");
+
+ hy_list.push_back("Ey_FEM_M180_0");
+ hy_list.push_back("Ey_FEM_M150_0");
+ hy_list.push_back("Ey_FEM_M120_0");
+ hy_list.push_back("Ey_FEM_M90_0");
+ hy_list.push_back("Ey_FEM_M60_0");
+ hy_list.push_back("Ey_FEM_M30_0");
+ hy_list.push_back("Ey_FEM_0_0");
+ hy_list.push_back("Ey_FEM_30_0");
+ hy_list.push_back("Ey_FEM_70_0");
+
+ int iVFD = 0;
+ float dVFD1 = 0;
+ float dVFD2 = 0;
+ for ( iVFD=0; iVFD<8; iVFD++ ){ // 2 of 9 will be seleted.
+ dVFD1 = m_VD - VFD0[iVFD];
+ dVFD2 = VFD0[iVFD+1] - m_VD;
+ if( dVFD2 >= 0 ) break;
+ }
+
+ h_Ex1 = static_cast<TH2F *>(hfile->Get((hx_list.at(iVFD)).c_str()));
+ h_Ex2 = static_cast<TH2F *>(hfile->Get((hx_list.at(iVFD+1)).c_str()));
+ h_Ey1 = static_cast<TH2F *>(hfile->Get((hy_list.at(iVFD)).c_str()));
+ h_Ey2 = static_cast<TH2F *>(hfile->Get((hy_list.at(iVFD+1)).c_str()));
+ h_ExHV = static_cast<TH2F *>(hfile->Get("Ex_FEM_0_100"));
+ h_EyHV = static_cast<TH2F *>(hfile->Get("Ey_FEM_0_100"));
+
+ float scalingFactor = m_VB / 100;
+
+ // For Electric field
+ for (int ix=0 ; ix <17 ; ix++){
+ for (int iy=0 ; iy<115 ; iy++) {
+ float Ex1 = h_Ex1->GetBinContent(ix+1,iy+1);
+ float Ex2 = h_Ex2->GetBinContent(ix+1,iy+1);
+ float ExHV = h_ExHV-> GetBinContent(ix+1,iy+1);
+ m_ExValue[ix][iy] = (Ex1*dVFD2+Ex2*dVFD1)/(dVFD1+dVFD2);
+ m_ExValue[ix][iy] += ExHV*scalingFactor;
+
+ float Ey1 = h_Ey1->GetBinContent(ix+1,iy+1);
+ float Ey2 = h_Ey2->GetBinContent(ix+1,iy+1);
+ float EyHV = h_EyHV-> GetBinContent(ix+1,iy+1);
+ m_EyValue[ix][iy] = (Ey1*dVFD2+Ey2*dVFD1)/(dVFD1+dVFD2);
+ m_EyValue[ix][iy] += EyHV*scalingFactor;
+ }
+ }
+
+ h_Ex1->Delete();
+ h_Ex2->Delete();
+ h_Ey1->Delete();
+ h_Ey2->Delete();
+ h_ExHV->Delete();
+ h_EyHV->Delete();
+ }
+
+ hfile -> Close();
+
+}
+
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.h b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..e0f505758cd8481152de61e0ec0879bb773aa94d
--- /dev/null
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.h
@@ -0,0 +1,123 @@
+#ifndef SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
+#define SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
+
+//-----------------------------------------------
+// 2020.8.12 Implementation in Athena by Manabu Togawa
+// 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)
+//-----------------------------------------------
+
+// C++ Standard Library
+#include <iostream>
+
+// ROOT
+#include "TH2.h"
+#include "TFile.h"
+
+// Athena
+#include "AthenaKernel/MsgStreamMember.h"
+#include "SiPropertiesSvc/ISiPropertiesSvc.h"
+#include "InDetConditionsSummaryService/ISiliconConditionsSvc.h"
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+#include "Identifier/IdentifierHash.h"
+#include "PathResolver/PathResolver.h"
+
+// Gaudi
+#include "GaudiKernel/PhysicalConstants.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// Random number
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGaussZiggurat.h" // for RandGaussZiggurat
+#include "CLHEP/Random/RandPoisson.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+namespace CLHEP {
+ class HepRandomEngine;
+}
+
+class SCT_InducedChargedModel {
+
+ public:
+
+ void Init(float vdepl, float vnbais, IdentifierHash m_hashId,
+ ServiceHandle<ISiliconConditionsSvc> m_siConditionsSvc);
+ void setRandomEngine(CLHEP::HepRandomEngine *rndmEngine) { m_rndmEngine = rndmEngine; };
+ void setMagneticField(double B){ m_B = - B*1000; }; // m_B in [Tesla]
+
+ void holeTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) ;
+ void electronTransport(double x0, double y0, double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2, IdentifierHash m_hashId, ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc ) ;
+
+ private:
+
+ void loadICMParameters();
+
+ bool electron(double x_e, double y_e, double &vx_e, double &vy_e, double &D_e,
+ IdentifierHash m_hashId,
+ ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) ;
+ bool hole(double x_h, double y_h, double &vx_h, double &vy_h, double &D_h,
+ IdentifierHash m_hashId,
+ ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc) ;
+ double induced (int istrip, double x, double y) ;
+ void EField( double x, double y, double &Ex, double &Ey ) ;
+
+
+ /// Log a message using the Athena controlled logging system
+ MsgStream& msg( MSG::Level lvl ) const { return m_msg << lvl; }
+ /// Check whether the logging system is active at the provided verbosity level
+ bool msgLvl( MSG::Level lvl ) { return m_msg.get().level() <= lvl; }
+
+ // Private message stream member
+ mutable Athena::MsgStreamMember m_msg;
+
+//------parameters given externally by jobOptions ------------------
+ int m_EfieldModel = 2 ; // 0(uniform E), 1(flat diode model), 2 (FEM solusions)
+ double m_VD = -30. ; // full depletion voltage [Volt] negative for type-P
+ double m_VB = 75. ; // applied bias voltage [Volt]
+ double m_T;
+ double m_B;
+ 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;
+
+//------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 = Gaudi::Units::k_Boltzmann / Gaudi::Units::joule; // [m^2*kg/s^2/K]
+ double m_e = Gaudi::Units::e_SI; // [Coulomb]
+
+//--------- parameters of induced charged model from root file ---------
+ TH2F *h_Potential_FDM;
+ TH2F *h_Potential_FEM;
+ TH2F *h_Ex1;
+ TH2F *h_Ex2;
+ TH2F *h_Ey1;
+ TH2F *h_Ey2;
+ TH2F *h_ExHV;
+ TH2F *h_EyHV;
+
+//---------- arrays of FEM analysis -----------------------------------
+// Storages for Ramo potential and Electric field.
+// In strip pitch directions :
+// Ramo potential : 80 divisions (81 points) with 5 um intervals from 40-440 um.
+// Electric field : 16 divisions (17 points) with 2.5 um intervals from 0-40 um.
+// In sensor depth directions (for both potential and Electric field):
+// 114 divisions (115 points) with 2.5 nm intervals for 285 um.
+
+ double m_PotentialValue[81][115];
+ double m_ExValue[17][115];
+ double m_EyValue[17][115];
+
+ CLHEP::HepRandomEngine *m_rndmEngine; //!< Random number generation engine
+
+};
+
+#endif // SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
index 9a1e828f694d2cda0633927718b4e3a8af485194..ced51d7d34ca85cf9d849f52caf375616b303f73 100755
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
@@ -63,6 +63,7 @@ parent)
m_doHistoTrap(false),
m_doRamo(false),
m_doCTrap(false),
+ m_doInducedChargedModel(false),
m_thistSvc(nullptr),
m_h_efieldz(nullptr),
m_h_efield(nullptr),
@@ -90,6 +91,7 @@ parent)
m_distortionsTool("SCT_DistortionsTool", this),
m_siConditionsSvc("SCT_SiliconConditionsSvc", name),
m_siPropertiesSvc("SCT_SiPropertiesSvc", name),
+ m_magFieldSvc("AtlasFieldSvc", name),
m_radDamageSvc("SCT_RadDamageSummarySvc", name),
m_element(0),
m_rndmEngine(0),
@@ -107,6 +109,7 @@ parent)
declareProperty("SmallStepLength", m_smallStepLength = 5);
declareProperty("SiConditionsSvc", m_siConditionsSvc);
declareProperty("SiPropertiesSvc", m_siPropertiesSvc);
+ declareProperty("MagFieldSvc", m_magFieldSvc);
// declareProperty("rndmEngineName",m_rndmEngineName="SCT_Digitization");
declareProperty("doDistortions", m_doDistortions,
"Simulation of module distortions");
@@ -124,6 +127,8 @@ parent)
"Tool to retrieve SCT distortions");
declareProperty("SCT_RadDamageSummarySvc", m_radDamageSvc);
declareProperty("isOverlay", m_isOverlay=false);
+ declareProperty("m_doInducedChargedModel", m_doInducedChargedModel,
+ "Flag for Induced Charged Model"); //
}
// Destructor:
@@ -152,6 +157,9 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
// Get ISCT_ModuleDistortionsTool
ATH_CHECK(m_distortionsTool.retrieve());
+ // Get IMagFieldSvc
+ ATH_CHECK(m_magFieldSvc.retrieve());
+
if (m_doTrapping) {
///////////////////////////////////////////////////
// -- Get Radiation Damage Service
@@ -236,6 +244,12 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
}
///////////////////////////////////////////////////
+ // Induced Charged Module. M.Togawa
+ if ( m_doInducedChargedModel ){
+ m_InducedChargedModel = std::make_unique<SCT_InducedChargedModel>();
+ m_InducedChargedModel->Init(m_vdepl,m_vbias,m_hashId,m_siConditionsSvc);
+ }
+
m_smallStepLength *= CLHEP::micrometer;
m_tSurfaceDrift *= CLHEP::ns;
@@ -295,7 +309,7 @@ float SCT_SurfaceChargesGenerator::DriftTime(float zhit) const {
float denominator = m_depletionVoltage + m_biasVoltage - (2.0 * zhit * m_depletionVoltage / m_thickness);
if (denominator <= 0.0) {
if (m_biasVoltage >= m_depletionVoltage) { // Should not happen
- if(!m_isOverlay) {
+ if (!m_isOverlay) {
ATH_MSG_ERROR("DriftTime: negative argument X for log(X) " << zhit);
}
return -1.0;
@@ -496,6 +510,21 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
float yhit = xPhi;
float zhit = xDep;
+ if ( m_doInducedChargedModel ){ // Setting magnetic field for the ICM.
+ HepGeom::Point3D<double> localpos(xhit,yhit,zhit);
+ HepGeom::Point3D<double> globalpos;
+ globalpos = m_element->globalPositionHit(localpos);
+ double point[3] = {globalpos.x(),globalpos.y(),globalpos.z()};
+ double field[3] = {0};
+ m_magFieldSvc->getField(point, field);
+ if ( m_isBarrel ){
+ m_InducedChargedModel->setMagneticField(field[2]);
+ }else{
+ double B_r = sqrt(pow(field[0],2)+pow(field[1],2));
+ m_InducedChargedModel->setMagneticField( B_r );
+ }
+ }
+
if (m_doDistortions) {
if (m_isBarrel) {// Only apply disortions to barrel
// modules
@@ -583,9 +612,12 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
if (tdrift > 0.0) {
float x1 = xhit + StepX * dstep;// (static_cast<float>(istep)+0.5) ;
float y1 = yhit + StepY * dstep;// (static_cast<float>(istep)+0.5) ;
- y1 += m_tanLorentz * zReadout; // !< Taking into account the magnetic
- // field
- float diffusionSigma = DiffusionSigma(zReadout);
+
+ float diffusionSigma = 0;
+ if ( !m_doInducedChargedModel ){
+ diffusionSigma = DiffusionSigma(zReadout);
+ y1 += m_tanLorentz * zReadout; // !< Taking into account the magnetic field
+ } // Should be treated in Induced Charged Model.
for (int i = 0; i < m_numberOfCharges; ++i) {
float rx = CLHEP::RandGaussZiggurat::shoot(m_rndmEngine);
@@ -733,6 +765,49 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
} // m_doTrapping==true
if (!m_doRamo) {
+
+ if ( m_doInducedChargedModel ){ // Induced Charged Model
+
+ // Charges storages for 50 ns. 0.5 ns steps.
+ double Q_m2[100]={0}, Q_m1[100]={0}, Q_00[100]={0}, Q_p1[100]={0},Q_p2[100]={0};
+
+ // Unit for y and z : mm -> cm in SCT_InducedChargedModel
+ m_InducedChargedModel->holeTransport(y0/10,z0/10,Q_m2,Q_m1,
+ Q_00,Q_p1,Q_p2,
+ m_hashId, m_siPropertiesSvc);
+ m_InducedChargedModel->electronTransport(y0/10,z0/10,Q_m2,
+ Q_m1,Q_00,Q_p1,Q_p2,
+ m_hashId, m_siPropertiesSvc);
+
+ for (int it=0; it<100; it++){
+ if ( Q_00[it] == 0.0) continue;
+
+ double ICM_time = (it+0.5)*0.5 + timeOfFlight;
+
+ double Q_new[5]={0};
+ Q_new[0] = Q_m2[it];
+ Q_new[1] = Q_m1[it];
+ Q_new[2] = Q_00[it];
+ Q_new[3] = Q_p1[it];
+ Q_new[4] = Q_p2[it];
+
+ for (int strip = -2; strip <= 2; strip++) {
+ double ystrip = y1 + strip * stripPitch;
+ SiLocalPosition position(m_element->hitLocalToLocal(x1, ystrip));
+ if (m_design->inActiveArea(position)) {
+ inserter(SiSurfaceCharge(position,
+ SiCharge(q1 * Q_new[strip+2],
+ ICM_time, hitproc,
+ trklink)));
+ }
+ }
+ }
+
+ } // End Induced Charged Model
+
+
+ if ( !m_doInducedChargedModel ){ // SCT_Digitization
+
SiLocalPosition position(m_element->hitLocalToLocal(xd,
yd));
if (m_design->inActiveArea(position)) {
@@ -772,6 +847,9 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
// "<<position<<" of the element is out of active area,
// charge = "<<q1) ;
}
+
+ } // End of SCT_Digitization
+
} // end of loop on charges
}
}
@@ -852,7 +930,7 @@ void SCT_SurfaceChargesGenerator::setVariables() {
m_electronHolePairsPerEnergy = m_siPropertiesSvc->getSiProperties(m_hashId).electronHolePairsPerEnergy();
// get sensor thickness and tg lorentz from SiDetectorDesign
- if(m_design) {
+ if (m_design) {
m_thickness = m_design->thickness();
} else {
m_thickness = 0.;
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
index 9dd7f02ba2a9d7f238ae6f174f13e6c164d64025..33f19ff4ce5fe3c72902165e0ae6113acd673f1e 100755
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
@@ -35,9 +35,13 @@
#include "SCT_Digitization/ISCT_SurfaceChargesGenerator.h"
#include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
+#include "SCT_InducedChargedModel.h"
+
#include "GaudiKernel/ToolHandle.h"
#include <iostream>
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+
// Charges and hits
class SiHit;
#include "Identifier/IdentifierHash.h"
@@ -87,7 +91,12 @@ private:
void setFixedTime(float fixedTime) {m_tfix = fixedTime;}
void setCosmicsRun(bool cosmicsRun) {m_cosmicsRun = cosmicsRun;}
void setComTimeFlag(bool useComTime) {m_useComTime = useComTime;}
- void setRandomEngine(CLHEP::HepRandomEngine *rndmEngine) {m_rndmEngine = rndmEngine;}
+ void setRandomEngine(CLHEP::HepRandomEngine *rndmEngine) {
+ m_rndmEngine = rndmEngine;
+ if ( m_doInducedChargedModel ){
+ m_InducedChargedModel->setRandomEngine( m_rndmEngine );
+ }
+ }
void setDetectorElement(const InDetDD::SiDetectorElement *ele) {m_element = ele; setVariables();}
/** create a list of surface charges from a hit */
@@ -162,6 +171,7 @@ private:
ServiceHandle<ISiliconConditionsSvc> m_siConditionsSvc;
ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc;
ServiceHandle<ISCT_RadDamageSummarySvc> m_radDamageSvc;
+ ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc;
const InDetDD::SiDetectorElement * m_element;
CLHEP::HepRandomEngine * m_rndmEngine; //!< Random Engine
@@ -180,6 +190,10 @@ private:
double m_center;
double m_tanLorentz;
bool m_isBarrel;
+
+ // For Induced Charged Module, M.Togawa
+ std::unique_ptr<SCT_InducedChargedModel> m_InducedChargedModel;
+ bool m_doInducedChargedModel; //!< Flag to use Induced Charged Model
};
#endif // SCT_SURFACECHARGESGENERATOR_H