diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt b/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
index cb3be644e67b7a359acaac4ef79d109ad0da9143..4209b461b0fde0d9a82ab895b31c4740ed35a4cf 100644
--- a/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
@@ -14,7 +14,7 @@ atlas_add_component( SCT_Digitization
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel CxxUtils PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesToolLib InDetIdentifier InDetReadoutGeometry SCT_ReadoutGeometry InDetSimData InDetConditionsSummaryService PathResolver SCT_ConditionsToolsLib SCT_ModuleDistortionsLib )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel CxxUtils PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesToolLib InDetIdentifier InDetReadoutGeometry SCT_ReadoutGeometry InDetSimData InDetConditionsSummaryService PathResolver SCT_ConditionsToolsLib SCT_ModuleDistortionsLib MagFieldConditions )
atlas_add_test( SCT_DigitizationMT_test
SCRIPT Digi_tf.py --inputHITSFile /cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/DigitizationTests/HITS.04919495._001041.pool.root.1 --conditionsTag default:OFLCOND-RUN12-SDR-25 --digiSeedOffset1 170 --digiSeedOffset2 170 --geometryVersion ATLAS-R2-2015-03-01-00 --DataRunNumber 222525 --outputRDOFile mc15_2015_ttbar.RDO.pool.root --preInclude HITtoRDO:SimulationJobOptions/preInclude.SCTOnlyConfig.py,Digitization/ForceUseOfAlgorithms.py --postInclude Digitization/FixDataDependenciesForMT.py --skipEvents 0 --maxEvents 100 --athenaopts=--threads=10
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.cxx b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4521adb2b3395deeaf4509b155d01ba43e11c5ce
--- /dev/null
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.cxx
@@ -0,0 +1,483 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+//-----------------------------------------------
+// 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_InducedChargeModel.h"
+
+// ROOT
+#include "TH2.h"
+#include "TFile.h"
+
+// C++ Standard Library
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+const double SCT_InducedChargeModel::s_kB = Gaudi::Units::k_Boltzmann / Gaudi::Units::joule; // [m^2*kg/s^2/K]
+const double SCT_InducedChargeModel::s_e = Gaudi::Units::e_SI; // [Coulomb]
+// Sorted depletion voltage values
+const std::vector<float> SCT_InducedChargeModel::s_VFD0 =
+ { -180., -150., -120., -90., -60., -30., 0., 30., 70.};
+const std::vector<std::string> SCT_InducedChargeModel::s_VFD0str =
+ {"M180", "M150", "M120", "M90", "M60", "M30", "0", "30", "70"};
+
+SCT_InducedChargeModel::SCT_InducedChargeModel(size_t maxHash, EFieldModel model) :
+ m_msg("SCT_InducedChargeModel"),
+ m_EFieldModel(model)
+{
+ loadICMParameters();
+ m_data.resize(maxHash);
+}
+
+SCT_InducedChargeModel::SCT_InducedChargeModelData*
+SCT_InducedChargeModel::setWaferData(const float vdepl,
+ const float vbias,
+ const InDetDD::SiDetectorElement* element,
+ const AtlasFieldCacheCondObj* fieldCondObj,
+ const ToolHandle<ISiliconConditionsTool> siConditionsTool,
+ CLHEP::HepRandomEngine* rndmEngine) const {
+ std::lock_guard<std::mutex> lock(m_mutex);
+
+ // If cache exists, cache is used.
+ SCT_InducedChargeModelData* p_data = m_data[element->identifyHash()].get();
+ if (p_data) return p_data;
+
+ ATH_MSG_DEBUG("--- Induced Charged Model Paramter (Begin) --------");
+
+ HepGeom::Point3D<double> localpos(0., 0., 0.); // The center of wafer
+ HepGeom::Point3D<double> globalpos = element->globalPositionHit(localpos);
+ double point[3] = {globalpos.x(), globalpos.y(), globalpos.z()};
+ double field[3] = {0., 0., 0.};
+ MagField::AtlasFieldCache fieldCache;
+ fieldCondObj->getInitializedCache(fieldCache);
+ fieldCache.getField(point, field);
+ const Amg::Vector3D magneticField(field[0], field[1], field[2]); // in kTesla
+
+ //---------------setting basic parameters---------------------------
+ std::unique_ptr<SCT_InducedChargeModelData> data =
+ std::make_unique<SCT_InducedChargeModelData>(vdepl,
+ vbias,
+ element,
+ magneticField,
+ m_bulk_depth,
+ m_EFieldModel,
+ siConditionsTool,
+ rndmEngine);
+
+ //--- set electric fields ---
+ setEField(*data);
+
+ //--------------------------------------------------------
+
+ ATH_MSG_DEBUG(" EFieldModel 0 (Flat Diode Model), 1 (FEM solusions) 2 (uniform E)\t"<< data->m_EFieldModel);
+ ATH_MSG_DEBUG(" DepletionVoltage_VD\t"<< data->m_VD <<" V");
+ ATH_MSG_DEBUG(" BiasVoltage_VB \t"<< data->m_VB <<" V");
+ ATH_MSG_DEBUG(" MagneticField_B \t"
+ << data->m_magneticField[Amg::x] << " "
+ << data->m_magneticField[Amg::y] << " "
+ << data->m_magneticField[Amg::z] <<" kTesla");
+ ATH_MSG_DEBUG(" Temperature_T \t"<< data->m_T <<" K");
+ ATH_MSG_DEBUG(" TransportTimeStep \t"<< m_transportTimeStep <<" ns");
+ ATH_MSG_DEBUG(" TransportTimeMax\t"<< m_transportTimeMax <<" ns");
+ ATH_MSG_DEBUG(" BulkDepth\t\t"<< m_bulk_depth << " cm");
+ ATH_MSG_DEBUG(" DepletionDepth\t\t"<< data->m_depletion_depth << " cm");
+ ATH_MSG_DEBUG(" StripPitch\t\t"<< m_strip_pitch << " cm");
+ ATH_MSG_DEBUG("--- Induced Charged Model Paramter (End) ---------");
+
+ m_data[element->identifyHash()] = std::move(data);
+
+ return m_data[element->identifyHash()].get();
+}
+
+//---------------------------------------------------------------------
+// transport
+//---------------------------------------------------------------------
+void SCT_InducedChargeModel::transport(const SCT_InducedChargeModelData& data,
+ const bool isElectron,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const {
+ // transport electrons and 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] NTransportSteps=100
+ // means 50 ns storages of charges in each strips.
+
+ double x = x0; // original electron/hole position [cm]
+ double y = y0; // original electron/hole position [cm]
+ bool isInBulk = true;
+ double t_current = 0.;
+ double qstrip[NStrips]; // induced charges on strips due to an electron or a hole
+ double vx, vy, D;
+ for (int istrip=StartStrip; istrip<=EndStrip; istrip++) {
+ qstrip[istrip+Offset] = (isElectron ? -1. : +1.) * induced(data, istrip, x, y); // original induced charge by electron or hole
+ }
+ while (t_current < m_transportTimeMax) {
+ if (!isInBulk) break;
+ if (!getVxVyD(data, isElectron, x, y, vx, vy, D, hashId, siPropertiesTool)) break;
+ double delta_y = vy * m_transportTimeStep / Gaudi::Units::second; // ns -> s
+ y += delta_y;
+ double dt = m_transportTimeStep; // [nsec]
+ if (isElectron) {
+ 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;
+ }
+ } else {
+ 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 += dt;
+ x += vx * dt / Gaudi::Units::second; // ns -> s
+ double diffusion = std::sqrt(2.* D * dt / Gaudi::Units::second); // ns -> s
+ y += diffusion * CLHEP::RandGaussZiggurat::shoot(data.m_rndmEngine, 0.0, 1.0);
+ x += diffusion * CLHEP::RandGaussZiggurat::shoot(data.m_rndmEngine, 0.0, 1.0);
+ if (isElectron) {
+ if (y < m_y_origin_min) {
+ y = m_y_origin_min;
+ isInBulk = false;
+ }
+ } else {
+ if (y > m_bulk_depth) {
+ y = m_bulk_depth;
+ isInBulk = false; // outside bulk region
+ }
+ }
+
+ // get induced current by subtracting induced charges
+ for (int istrip=StartStrip; istrip<=EndStrip; istrip++) {
+ double qnew = (isElectron ? -1. : +1.) * induced(data, istrip, x, y);
+ int jj = istrip + Offset;
+ double dq = qnew - qstrip[jj];
+ qstrip[jj] = qnew;
+
+ int jt = static_cast<int>((t_current+0.001) / m_transportTimeStep);
+ if (jt < NTransportSteps) {
+ 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 or hole tracing
+}
+
+//---------------------------------------------------------------------
+// holeTransport
+//---------------------------------------------------------------------
+void SCT_InducedChargeModel::holeTransport(const SCT_InducedChargeModelData& data,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const {
+ // Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100] NTransportSteps=100
+ const bool isElectron = false;
+ transport(data,
+ isElectron,
+ x0, y0,
+ Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+ hashId,
+ siPropertiesTool);
+}
+
+//---------------------------------------------------------------------
+// electronTransport
+//---------------------------------------------------------------------
+void SCT_InducedChargeModel::electronTransport(const SCT_InducedChargeModelData& data,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const {
+ // Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100] NTransportSteps=100
+ const bool isElectron = true;
+ transport(data,
+ isElectron,
+ x0, y0,
+ Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+ hashId,
+ siPropertiesTool);
+}
+
+//---------------------------------------------------------------
+// Get parameters for electron or hole transport
+//---------------------------------------------------------------
+bool SCT_InducedChargeModel::getVxVyD(const SCT_InducedChargeModelData& data,
+ const bool isElectron,
+ const double x, const double y, double& vx, double& vy, double& D,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const {
+ double Ex, Ey;
+ getEField(data, x, y, Ex, Ey); // [V/cm]
+ if (Ey > 0.) {
+ if (isElectron) {
+ // because electron has negative charge
+ Ex *= -1.;
+ Ey *= -1.;
+ }
+ const double E = std::sqrt(Ex*Ex+Ey*Ey);
+ const double mu = (isElectron ?
+ siPropertiesTool->getSiProperties(hashId).calcElectronDriftMobility(data.m_T, E*CLHEP::volt/CLHEP::cm) :
+ siPropertiesTool->getSiProperties(hashId).calcHoleDriftMobility (data.m_T, E*CLHEP::volt/CLHEP::cm))
+ * (CLHEP::volt/CLHEP::cm)/(CLHEP::cm/CLHEP::s);
+ const double v = mu * E;
+ const double r = (isElectron ?
+ siPropertiesTool->getSiProperties(hashId).calcElectronHallFactor(data.m_T) :
+ siPropertiesTool->getSiProperties(hashId).calcHoleHallFactor(data.m_T));
+
+ const double tanLA = data.m_element->design().readoutSide()
+ * r * mu * (isElectron ? -1. : +1.) // because e has negative charge.
+ * data.m_element->hitDepthDirection()
+ * data.m_element->hitPhiDirection()
+ * (data.m_element->normal().cross(data.m_magneticField)).dot(data.m_element->phiAxis())
+ / Gaudi::Units::tesla * Gaudi::Units::gauss * 1000.; // 1000. is to convert kTesla to Tesla.
+
+ const double secLA = std::sqrt(1.+tanLA*tanLA);
+ const double cosLA = 1./secLA;
+ const double sinLA = tanLA/secLA;
+ vy = v * (Ey*cosLA - Ex*sinLA)/E;
+ vx = v * (Ex*cosLA + Ey*sinLA)/E;
+ D = s_kB * data.m_T * mu/ s_e;
+ return true;
+ }
+ return false;
+}
+
+//-------------------------------------------------------------------
+// calculation of induced charge using Weighting (Ramo) function
+//-------------------------------------------------------------------
+double SCT_InducedChargeModel::induced(const SCT_InducedChargeModelData& data,
+ const int istrip, const double x, const double y) const
+{
+ // 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]
+
+ // Ramo potential : 80 divisions (81 points) with 5 um intervals from 40-440 um.
+ // In sensor depth directions 114 divisions (115 points) with 2.5 nm intervals for 285 um.
+ const 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 = std::abs(x-xc);
+ int ix = static_cast<int>(dx / deltax);
+ if (ix >= NRamoPoints) return 0.;
+ int iy = static_cast<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[data.m_EFieldModel][ix ][iy ] *(1.-fx)*(1.-fy)
+ + m_PotentialValue[data.m_EFieldModel][ix1][iy ] * fx *(1.-fy)
+ + m_PotentialValue[data.m_EFieldModel][ix ][iy1] *(1.-fx)* fy
+ + m_PotentialValue[data.m_EFieldModel][ix1][iy1] * fx* fy;
+
+ return P;
+}
+
+//--------------------------------------------------------------
+// Electric Field Ex and Ey
+//--------------------------------------------------------------
+void SCT_InducedChargeModel::getEField(const SCT_InducedChargeModelData& data,
+ const double x, const double y, double& Ex, double& Ey) const {
+ // m_EFieldModel == FlatDiodeModel 0; flat diode model
+ // m_EFieldModel == FEMsolutions 1; FEM solusions
+ // m_EFieldModel == UniformE 2; uniform E-field model
+ // x == 0.0040 [cm] : at the center of strip
+ // x must be within 0 and 0.008 [cm]
+
+ // Electric field : 16 divisions (17 points) with 2.5 um intervals from 0-40 um.
+ // In sensor depth directions 114 divisions (115 points) with 2.5 nm intervals for 285 um.
+ const 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 (data.m_EFieldModel==FEMsolutions) {
+ int iy = static_cast<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 = static_cast<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 = data.m_ExValue[ix][iy]; Ex10 = data.m_ExValue[ix1][iy];
+ Ex01 = data.m_ExValue[ix][iy1]; Ex11 = data.m_ExValue[ix1][iy1];
+ Ey00 = data.m_EyValue[ix][iy]; Ey10 = data.m_EyValue[ix1][iy];
+ Ey01 = data.m_EyValue[ix][iy1]; Ey11 = data.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;
+ }
+ //---------- case for uniform electric field ------------------------
+ else if (data.m_EFieldModel==UniformE) {
+ if (m_bulk_depth - y < data.m_depletion_depth && data.m_depletion_depth >0.) {
+ Ey = data.m_VB / data.m_depletion_depth;
+ } else {
+ Ey = 0.;
+ }
+ }
+ //---------- case for flat diode model ------------------------------
+ else if (data.m_EFieldModel==FlatDiodeModel) {
+ if (data.m_VB > std::abs(data.m_VD)) {
+ Ey = (data.m_VB+data.m_VD)/m_bulk_depth - 2.*data.m_VD*(m_bulk_depth-y)/(m_bulk_depth*m_bulk_depth);
+ } else {
+ if (m_bulk_depth - y < data.m_depletion_depth && data.m_depletion_depth >0.) {
+ double Emax = 2.* data.m_depletion_depth * std::abs(data.m_VD) / (m_bulk_depth*m_bulk_depth);
+ Ey = Emax*(1.-(m_bulk_depth-y)/data.m_depletion_depth);
+ } else {
+ Ey = 0.;
+ }
+ }
+ }
+}
+
+/////////////////////////////////////////////////////////////////////
+size_t SCT_InducedChargeModel::getFEMIndex(SCT_InducedChargeModelData& data) const {
+ // Return index for s_VFD0 and EFieldModel.
+ // If vdepl is out of range of s_VFD0, EFieldModel of FlatDiodeModel is returned.
+ size_t iVFD = 0;
+ if (data.m_VD<s_VFD0[0] || data.m_VD>s_VFD0[s_VFD0.size()-1]) { // Out of range
+ data.m_EFieldModel = FlatDiodeModel;
+ } else {
+ float dVFD2 = 0.;
+ for (iVFD=0; iVFD<s_VFD0.size()-1; iVFD++) {
+ dVFD2 = s_VFD0[iVFD+1] - data.m_VD;
+ if (dVFD2 >= 0.) break;
+ }
+ }
+ return iVFD;
+}
+
+/////////////////////////////////////////////////////////////////////
+
+void SCT_InducedChargeModel::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());
+
+ // For Ramo Potential
+ TH2F* h_Potential_FDM = static_cast<TH2F*>(hfile->Get("Potential_FDM"));
+ TH2F* h_Potential_FEM = static_cast<TH2F*>(hfile->Get("Potential_FEM"));
+ // FDM case keeps only FDM potential.
+ // FEM case keeps both FDM and FEM potentials.
+ if (m_EFieldModel==FlatDiodeModel || m_EFieldModel==FEMsolutions) {
+ if (m_EFieldModel==FlatDiodeModel) m_PotentialValue.resize(FlatDiodeModel+1);
+ else m_PotentialValue.resize(FEMsolutions+1);
+ for (int ix=0; ix<NRamoPoints; ix++) {
+ for (int iy=0; iy<NDepthPoints; iy++) {
+ m_PotentialValue[FlatDiodeModel][ix][iy] = h_Potential_FDM->GetBinContent(ix+1, iy+1);
+ if (m_EFieldModel==FEMsolutions) {
+ m_PotentialValue[FEMsolutions][ix][iy] = h_Potential_FEM->GetBinContent(ix+1, iy+1);
+ }
+ }
+ }
+ }
+ h_Potential_FDM->Delete();
+ h_Potential_FEM->Delete();
+
+ if (m_EFieldModel==FEMsolutions) {
+ size_t i=0;
+ m_ExValue.resize(s_VFD0str.size()+1);
+ m_EyValue.resize(s_VFD0str.size()+1);
+ for (const std::string& str : s_VFD0str) {
+ TH2F* h_Ex = static_cast<TH2F*>(hfile->Get(("Ex_FEM_"+str+"_0").c_str()));
+ TH2F* h_Ey = static_cast<TH2F*>(hfile->Get(("Ey_FEM_"+str+"_0").c_str()));
+ for (int ix=0; ix <NEFieldPoints; ix++){
+ for (int iy=0; iy<NDepthPoints; iy++) {
+ m_ExValue[i][ix][iy] = h_Ex->GetBinContent(ix+1, iy+1);
+ m_EyValue[i][ix][iy] = h_Ey->GetBinContent(ix+1, iy+1);
+ }
+ }
+ h_Ex->Delete();
+ h_Ey->Delete();
+ i++;
+ }
+
+ TH2F* h_ExHV = static_cast<TH2F*>(hfile->Get("Ex_FEM_0_100"));
+ TH2F* h_EyHV = static_cast<TH2F*>(hfile->Get("Ey_FEM_0_100"));
+ for (int ix=0; ix <NEFieldPoints; ix++){
+ for (int iy=0; iy<NDepthPoints; iy++) {
+ m_ExValue[i][ix][iy] = h_ExHV->GetBinContent(ix+1, iy+1);
+ m_EyValue[i][ix][iy] = h_EyHV->GetBinContent(ix+1, iy+1);
+ }
+ }
+ h_ExHV->Delete();
+ h_EyHV->Delete();
+ }
+
+ hfile->Close();
+}
+
+void SCT_InducedChargeModel::setEField(SCT_InducedChargeModelData& data) const {
+ if (data.m_EFieldModel!=FEMsolutions) return;
+
+ size_t iVFD = getFEMIndex(data);
+ if (data.m_EFieldModel==FlatDiodeModel) {
+ ATH_MSG_INFO("Changed to Flat Diode Model since deplettion volage is out of range. ("
+ << s_VFD0[0] << " < m_VD < " << s_VFD0[s_VFD0.size()-1] << " is allowed.)");
+ return;
+ }
+
+ const float dVFD1 = data.m_VD - s_VFD0[iVFD];
+ const float dVFD2 = s_VFD0[iVFD+1] - data.m_VD;
+
+ const float scalingFactor = data.m_VB / 100.; // Reference voltage is 100 V
+
+ const size_t iHV = s_VFD0.size();
+
+ // For Electric field
+ for (int ix=0; ix <NEFieldPoints; ix++){
+ for (int iy=0; iy<NDepthPoints; iy++) {
+ float Ex1 = m_ExValue[iVFD ][ix][iy];
+ float Ex2 = m_ExValue[iVFD+1][ix][iy];
+ float ExHV = m_ExValue[iHV ][ix][iy];
+ data.m_ExValue[ix][iy] = (Ex1*dVFD2+Ex2*dVFD1)/(dVFD1+dVFD2);
+ data.m_ExValue[ix][iy] += ExHV*scalingFactor;
+
+ float Ey1 = m_EyValue[iVFD ][ix][iy];
+ float Ey2 = m_EyValue[iVFD+1][ix][iy];
+ float EyHV = m_EyValue[iHV ][ix][iy];
+ data.m_EyValue[ix][iy] = (Ey1*dVFD2+Ey2*dVFD1)/(dVFD1+dVFD2);
+ data.m_EyValue[ix][iy] += EyHV*scalingFactor;
+ }
+ }
+}
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.h b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..d25c00fa6d2c26c9f9f3a2979bd49c1a0a70b36a
--- /dev/null
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.h
@@ -0,0 +1,183 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#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)
+//-----------------------------------------------
+
+// Athena
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "AthenaKernel/MsgStreamMember.h"
+#include "Identifier/IdentifierHash.h"
+#include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
+#include "InDetReadoutGeometry/SiDetectorElement.h"
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
+#include "PathResolver/PathResolver.h"
+#include "SiPropertiesTool/ISiPropertiesTool.h"
+
+// Gaudi
+#include "GaudiKernel/PhysicalConstants.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// Random number
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGaussZiggurat.h" // for RandGaussZiggurat
+#include "CLHEP/Random/RandPoisson.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+// C++ Standard Library
+#include <array>
+#include <memory>
+#include <mutex>
+#include <string>
+#include <utility>
+#include <vector>
+
+class SCT_InducedChargeModel {
+
+ public:
+ enum EFieldModel {FlatDiodeModel=0, FEMsolutions=1, UniformE=2};
+ enum TransportStep {NTransportSteps=100};
+ enum FEMNums {NRamoPoints=81, NEFieldPoints=17, NDepthPoints=115};
+ enum InducedStrips {StartStrip=-2, EndStrip=+2, Offset=-StartStrip, NStrips=EndStrip+Offset+1};
+
+ struct SCT_InducedChargeModelData {
+ float m_VD; // full depletion voltage [Volt] negative for type-P
+ float m_VB; // applied bias voltage [Volt]
+ float m_T; // temperature
+ float m_depletion_depth;
+ const InDetDD::SiDetectorElement* m_element;
+ Amg::Vector3D m_magneticField;
+ EFieldModel m_EFieldModel;
+ CLHEP::HepRandomEngine* m_rndmEngine;
+ std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_ExValue;
+ std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_EyValue;
+
+ SCT_InducedChargeModelData(const float vdepl,
+ const float vbias,
+ const InDetDD::SiDetectorElement* element,
+ const Amg::Vector3D& magneticField, // in kTesla
+ const float bulk_depth,
+ const EFieldModel model,
+ const ToolHandle<ISiliconConditionsTool> siConditionsTool,
+ CLHEP::HepRandomEngine* rndmEngine) {
+ m_VD = vdepl; // full depletion voltage [Volt] negative for type-P
+ m_VB = vbias; // applied bias voltage [Volt]
+ m_element = element;
+ m_magneticField = magneticField;
+
+ //------------ find delepletion deph for model=0 and 1 -------------
+ m_depletion_depth = bulk_depth;
+ // for type N (before type inversion)
+ if (m_VD >= 0.) {
+ if (m_VB < m_VD) m_depletion_depth = std::sqrt(m_VB/m_VD) * bulk_depth;
+ } else {
+ // for type P
+ if (m_VB <= std::abs(m_VD)) m_depletion_depth = 0.;
+ }
+
+ m_EFieldModel = model;
+ m_T = siConditionsTool->temperature(m_element->identifyHash()) + Gaudi::Units::STP_Temperature;
+ m_rndmEngine = rndmEngine;
+ }
+ };
+
+ SCT_InducedChargeModel(size_t maxHash, EFieldModel model=FEMsolutions);
+
+ SCT_InducedChargeModelData*
+ setWaferData(const float vdepl,
+ const float vbias,
+ const InDetDD::SiDetectorElement* element,
+ const AtlasFieldCacheCondObj* fieldCondObj,
+ const ToolHandle<ISiliconConditionsTool> siConditionsTool,
+ CLHEP::HepRandomEngine* rndmEngine) const;
+
+ void setEField(SCT_InducedChargeModelData& data) const;
+
+ void transport(const SCT_InducedChargeModelData& data,
+ const bool isElectron,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+ void holeTransport(const SCT_InducedChargeModelData& data,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+ void electronTransport(const SCT_InducedChargeModelData& data,
+ const double x0, const double y0,
+ double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+
+ private:
+
+ void loadICMParameters();
+
+ bool getVxVyD(const SCT_InducedChargeModelData& data,
+ const bool isElectron,
+ const double x, const double y, double& vx, double& vy, double& D,
+ const IdentifierHash hashId,
+ const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+ double induced(const SCT_InducedChargeModelData& data,
+ const int istrip, const double x, const double y) const;
+ void getEField(const SCT_InducedChargeModelData& data,
+ const double x, const double y, double& Ex, double& Ey) const;
+
+ size_t getFEMIndex(SCT_InducedChargeModelData& data) const;
+
+ /// 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) const { return m_msg.get().level() <= lvl; }
+
+ //-------- parameters for e, h transport --------------------------------
+ static const double s_kB; // [m^2*kg/s^2/K]
+ static const double s_e; // [Coulomb]
+
+ // Private message stream member
+ mutable Athena::MsgStreamMember m_msg ATLAS_THREAD_SAFE;
+
+ //------parameters given externally by jobOptions ------------------
+ EFieldModel m_EFieldModel; // 0 (flat diode model), 1 (FEM solusions), 2 (uniform E)
+ double m_transportTimeStep = 0.50; // one step side in time [nsec]
+ double m_transportTimeMax = m_transportTimeStep*NTransportSteps; // maximun tracing time [nsec]
+
+ //------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_y_origin_min = 0.;
+
+ //---------- 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.
+
+ std::vector<std::array<std::array<double, NDepthPoints>, NRamoPoints>> m_PotentialValue;
+ std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_ExValue;
+ std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_EyValue;
+
+ // Cache of SCT_InducedChargeModelData for each wafer.
+ // Assuming wafer parameters do not change during a job.
+ // This should be almost always satisfied.
+ mutable std::vector<std::unique_ptr<SCT_InducedChargeModelData>> m_data ATLAS_THREAD_SAFE; // Guarded by m_mutex
+ // To protect concurrent access to m_data
+ mutable std::mutex m_mutex;
+
+ static const std::vector<float> s_VFD0;
+ static const std::vector<std::string> s_VFD0str;
+};
+
+#endif // SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx
deleted file mode 100644
index 3f326f0814a2586fc062a8718f273a33802479cf..0000000000000000000000000000000000000000
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.cxx
+++ /dev/null
@@ -1,410 +0,0 @@
-/*
- Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
-*/
-
-//-----------------------------------------------
-// 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"
-
-// ROOT
-#include "TH2.h"
-#include "TFile.h"
-
-// C++ Standard Library
-#include <algorithm>
-#include <cmath>
-#include <iostream>
-
-const double SCT_InducedChargedModel::s_kB = Gaudi::Units::k_Boltzmann / Gaudi::Units::joule; // [m^2*kg/s^2/K]
-const double SCT_InducedChargedModel::s_e = Gaudi::Units::e_SI; // [Coulomb]
-const std::vector<float> SCT_InducedChargedModel::s_VFD0 =
- { -180., -150., -120., -90., -60., -30., 0., 30., 70.};
-const std::vector<std::string> SCT_InducedChargedModel::s_VFD0str =
- {"M180", "M150", "M120", "M90", "M60", "M30", "0", "30", "70"};
-
-void SCT_InducedChargedModel::init(const float vdepl, const float vbias, const IdentifierHash hashId, ToolHandle<ISiliconConditionsTool>& siConditionsTool) {
-
- 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 = siConditionsTool->temperature(hashId) + Gaudi::Units::STP_Temperature;
-
- //------------------------ initialize subfunctions ------------
-
- loadICMParameters();
-
- //------------ find delepletion deph for model=0 and 1 -------------
- m_depletion_depth = m_bulk_depth;
- // for type N (before type inversion)
- if (m_VD >= 0.) {
- if (m_VB < m_VD) m_depletion_depth = std::sqrt(m_VB/m_VD) * m_bulk_depth;
- } else {
- // for type P
- if (m_VB <= std::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;
-}
-
-//---------------------------------------------------------------------
-// transport
-//---------------------------------------------------------------------
-void SCT_InducedChargedModel::transport(const bool isElectron,
- const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const {
- // transport electrons and 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] NTransportSteps=100
- // means 50 ns storages of charges in each strips.
-
- double x = x0; // original electron/hole position [cm]
- double y = y0; // original electron/hole position [cm]
- bool isInBulk = true;
- double t_current = 0.;
- double qstrip[NStrips]; // induced charges on strips due to an electron or a hole
- double vx, vy, D;
- for (int istrip=StartStrip; istrip<=EndStrip; istrip++) {
- qstrip[istrip+Offset] = (isElectron ? -1. : +1.) * induced(istrip, x, y); // original induced charge by electron or hole
- }
- while (t_current < m_transportTimeMax) {
- if (!isInBulk) break;
- if (!getVxVyD(isElectron, x, y, vx, vy, D, hashId, siPropertiesTool)) break;
- double delta_y = vy * m_transportTimeStep / Gaudi::Units::second; // ns -> s
- y += delta_y;
- double dt = m_transportTimeStep; // [nsec]
- if (isElectron) {
- 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;
- }
- } else {
- 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 += dt;
- x += vx * dt / Gaudi::Units::second; // ns -> s
- double diffusion = std::sqrt(2.* D * dt / Gaudi::Units::second); // ns -> s
- y += diffusion * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
- x += diffusion * CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, 0.0, 1.0);
- if (isElectron) {
- if (y < m_y_origin_min) {
- y = m_y_origin_min;
- isInBulk = false;
- }
- } else {
- if (y > m_bulk_depth) {
- y = m_bulk_depth;
- isInBulk = false; // outside bulk region
- }
- }
-
- // get induced current by subtracting induced charges
- for (int istrip=StartStrip; istrip<=EndStrip; istrip++) {
- double qnew = (isElectron ? -1. : +1.) * induced(istrip, x, y);
- int jj = istrip + Offset;
- double dq = qnew - qstrip[jj];
- qstrip[jj] = qnew;
-
- int jt = static_cast<int>((t_current+0.001) / m_transportTimeStep);
- if (jt < NTransportSteps) {
- 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 or hole tracing
-}
-
-//---------------------------------------------------------------------
-// holeTransport
-//---------------------------------------------------------------------
-void SCT_InducedChargedModel::holeTransport(const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const {
- // Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100] NTransportSteps=100
- const bool isElectron = false;
- transport(isElectron,
- x0, y0,
- Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
- hashId,
- siPropertiesTool);
-}
-
-//---------------------------------------------------------------------
-// electronTransport
-//---------------------------------------------------------------------
-void SCT_InducedChargedModel::electronTransport(const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const {
- // Q_m2[100],Q_m1[100],Q_00[100],Q_p1[100],Q_p2[100] NTransportSteps=100
- const bool isElectron = true;
- transport(isElectron,
- x0, y0,
- Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
- hashId,
- siPropertiesTool);
-}
-
-//---------------------------------------------------------------
-// Get parameters for electron or hole transport
-//---------------------------------------------------------------
-bool SCT_InducedChargedModel::getVxVyD(const bool isElectron,
- const double x, const double y, double& vx, double& vy, double& D,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const {
- double Ex, Ey;
- EField(x, y, Ex, Ey); // [V/cm]
- if (Ey > 0.) {
- if (isElectron) {
- // because electron has negative charge
- Ex *= -1.;
- Ey *= -1.;
- }
- const double E = std::sqrt(Ex*Ex+Ey*Ey);
- const double mu = (isElectron ?
- siPropertiesTool->getSiProperties(hashId).calcElectronDriftMobility(m_T, E*CLHEP::volt/CLHEP::cm) :
- siPropertiesTool->getSiProperties(hashId).calcHoleDriftMobility (m_T, E*CLHEP::volt/CLHEP::cm))
- * (CLHEP::volt/CLHEP::cm)/(CLHEP::cm/CLHEP::s);
- const double v = mu * E;
- const double r = (isElectron ?
- 1.13 + 0.0008*(m_T - Gaudi::Units::STP_Temperature) :
- 0.72 - 0.0005*(m_T - Gaudi::Units::STP_Temperature));
- const double tanLA = r * mu * (isElectron ? -1. : +1.) * m_B / Gaudi::Units::tesla * Gaudi::Units::gauss; // because e has negative charge.
- const double secLA = std::sqrt(1.+tanLA*tanLA);
- const double cosLA = 1./secLA;
- const double sinLA = tanLA/secLA;
- vy = v * (Ey*cosLA - Ex*sinLA)/E;
- vx = v * (Ex*cosLA + Ey*sinLA)/E;
- D = s_kB * m_T * mu/ s_e;
- return true;
- }
- return false;
-}
-
-//-------------------------------------------------------------------
-// calculation of induced charge using Weighting (Ramo) function
-//-------------------------------------------------------------------
-double SCT_InducedChargedModel::induced(const int istrip, const double x, const double y) const
-{
- // 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 = std::abs(x-xc);
- int ix = static_cast<int>(dx / deltax);
- if (ix >= NRamoPoints) return 0.;
- int iy = static_cast<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(const double x, const double y, double& Ex, double& Ey) const {
- // m_EFieldModel == UniformE 0; uniform E-field model
- // m_EFieldModel == FlatDiodeModel 1; flat diode model
- // m_EFieldModel == FEMsolutions 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==FEMsolutions) {
- int iy = static_cast<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 = static_cast<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;
- }
- //---------- case for uniform electric field ------------------------
- else if (m_EFieldModel==UniformE) {
- if (m_bulk_depth - y < m_depletion_depth && m_depletion_depth >0.) {
- Ey = m_VB / m_depletion_depth;
- } else {
- Ey = 0.;
- }
- }
- //---------- case for flat diode model ------------------------------
- else if (m_EFieldModel==FlatDiodeModel) {
- if (m_VB > std::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 * std::abs(m_VD) / (m_bulk_depth*m_bulk_depth);
- Ey = Emax*(1.-(m_bulk_depth-y)/m_depletion_depth);
- } else {
- Ey = 0.;
- }
- }
- }
-}
-
-/////////////////////////////////////////////////////////////////////
-
-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());
-
- TH2F* h_Potential_FDM = static_cast<TH2F*>(hfile->Get("Potential_FDM"));
- TH2F* h_Potential_FEM = static_cast<TH2F*>(hfile->Get("Potential_FEM"));
-
- const float minBiasV_ICM_FEM = *std::min_element(s_VFD0.begin(), s_VFD0.end());
- const float maxBiasV_ICM_FEM = *std::max_element(s_VFD0.begin(), s_VFD0.end());
- if (m_VD<minBiasV_ICM_FEM || m_VD>maxBiasV_ICM_FEM) {
- m_EFieldModel = FlatDiodeModel; // Change to FDM
- ATH_MSG_INFO("Changed to Flat Diode Model since deplettion volage is out of range. ("
- << minBiasV_ICM_FEM << " < m_VD < " << maxBiasV_ICM_FEM << " is allowed.)");
- }
-
- // For Ramo Potential
- for (int ix=0; ix<NRamoPoints; ix++) {
- for (int iy=0; iy<NDepthPoints; iy++) {
- if (m_EFieldModel==FEMsolutions) { // 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==FEMsolutions) {
-
- std::vector<std::string> hx_list;
- std::vector<std::string> hy_list;
- for (const std::string& str : s_VFD0str) {
- hx_list.push_back("Ex_FEM_"+str+"_0");
- hy_list.push_back("Ey_FEM_"+str+"_0");
- }
-
- size_t iVFD = 0;
- float dVFD1 = 0.;
- float dVFD2 = 0.;
- for (iVFD=0; iVFD<s_VFD0.size()-1; iVFD++) { // 2 of 9 will be seleted.
- dVFD1 = m_VD - s_VFD0[iVFD];
- dVFD2 = s_VFD0[iVFD+1] - m_VD;
- if (dVFD2 >= 0.) break;
- }
-
- TH2F* h_Ex1 = static_cast<TH2F*>(hfile->Get((hx_list.at(iVFD)).c_str()));
- TH2F* h_Ex2 = static_cast<TH2F*>(hfile->Get((hx_list.at(iVFD+1)).c_str()));
- TH2F* h_Ey1 = static_cast<TH2F*>(hfile->Get((hy_list.at(iVFD)).c_str()));
- TH2F* h_Ey2 = static_cast<TH2F*>(hfile->Get((hy_list.at(iVFD+1)).c_str()));
- TH2F* h_ExHV = static_cast<TH2F*>(hfile->Get("Ex_FEM_0_100"));
- TH2F* h_EyHV = static_cast<TH2F*>(hfile->Get("Ey_FEM_0_100"));
-
- float scalingFactor = m_VB / 100.; // Reference voltage is 100 V??
-
- // For Electric field
- for (int ix=0; ix <NEFieldPoints; ix++){
- for (int iy=0; iy<NDepthPoints; 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
deleted file mode 100644
index d302cd03c840e66249629c309e258dfa9109609b..0000000000000000000000000000000000000000
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.h
+++ /dev/null
@@ -1,123 +0,0 @@
-/*
- Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
-*/
-
-#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)
-//-----------------------------------------------
-
-// Athena
-#include "AthenaKernel/IAtRndmGenSvc.h"
-#include "AthenaKernel/MsgStreamMember.h"
-#include "CxxUtils/checker_macros.h"
-#include "Identifier/IdentifierHash.h"
-#include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
-#include "PathResolver/PathResolver.h"
-#include "SiPropertiesTool/ISiPropertiesTool.h"
-
-// Gaudi
-#include "GaudiKernel/PhysicalConstants.h"
-#include "GaudiKernel/ToolHandle.h"
-
-// Random number
-#include "CLHEP/Random/RandFlat.h"
-#include "CLHEP/Random/RandGaussZiggurat.h" // for RandGaussZiggurat
-#include "CLHEP/Random/RandPoisson.h"
-#include "CLHEP/Units/SystemOfUnits.h"
-
-// C++ Standard Library
-#include <string>
-#include <vector>
-
-class SCT_InducedChargedModel {
-
- public:
- enum EFieldModel {UniformE=0, FlatDiodeModel=1, FEMsolutions=2};
- enum TransportStep {NTransportSteps=100};
-
- void init(const float vdepl, const float vbias, const IdentifierHash hashId,
- ToolHandle<ISiliconConditionsTool>& siConditionsTool);
- void setRandomEngine(CLHEP::HepRandomEngine* rndmEngine) { m_rndmEngine = rndmEngine; };
- void setMagneticField(const double B) { m_B = - B*1000.; }; // m_B in [Tesla], B in kTesla
-
- void transport(const bool isElectron,
- const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
- void holeTransport(const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
- void electronTransport(const double x0, const double y0,
- double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
-
- private:
- enum InducedStrips {StartStrip=-2, EndStrip=+2, Offset=-StartStrip, NStrips=EndStrip+Offset+1};
-
- void loadICMParameters();
-
- bool getVxVyD(const bool isElectron,
- const double x, const double y, double& vx, double& vy, double& D,
- const IdentifierHash hashId,
- ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
- double induced(const int istrip, const double x, const double y) const;
- void EField(const double x, const double y, double& Ex, double& Ey) const;
-
-
- /// 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; }
-
- //-------- parameters for e, h transport --------------------------------
- static const double s_kB; // [m^2*kg/s^2/K]
- static const double s_e; // [Coulomb]
-
- // Private message stream member
- mutable Athena::MsgStreamMember m_msg ATLAS_THREAD_SAFE;
-
- CLHEP::HepRandomEngine* m_rndmEngine; //!< Random number generation engine
-
- //------parameters given externally by jobOptions ------------------
- EFieldModel m_EFieldModel = FEMsolutions; // 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 = m_transportTimeStep*NTransportSteps; // maximun tracing time [nsec]
-
- //------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 = 0.;
-
- //---------- 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.
-
- enum FEMNums {NRamoPoints=81, NEFieldPoints=17, NDepthPoints=115};
-
- double m_PotentialValue[NRamoPoints][NDepthPoints];
- double m_ExValue[NEFieldPoints][NDepthPoints];
- double m_EyValue[NEFieldPoints][NDepthPoints];
-
- static const std::vector<float> s_VFD0;
- static const std::vector<std::string> s_VFD0str;
-};
-
-#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 afbcb52c9653b4a8f36c67b5e80a348b7b717742..5e360577cedf7d48bcd86b7ab4ee1f4784a33428 100644
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
@@ -10,6 +10,7 @@
// Athena
#include "GeneratorObjects/HepMcParticleLink.h"
+#include "InDetIdentifier/SCT_ID.h"
#include "InDetSimEvent/SiHit.h" // for SiHit, SiHit::::xDep, etc
#include "HitManagement/TimedHitPtr.h" // for TimedHitPtr
@@ -49,6 +50,8 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
// Get ISiliconConditionsSvc
ATH_CHECK(m_siConditionsTool.retrieve());
+ ATH_CHECK(m_fieldCacheCondObjInputKey.initialize(m_doInducedChargeModel));
+
if (m_doTrapping) {
// -- Get Radiation Damage Tool
ATH_CHECK(m_radDamageTool.retrieve());
@@ -129,6 +132,13 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
}
///////////////////////////////////////////////////
+ // Induced Charge Module.
+ if (m_doInducedChargeModel) {
+ const SCT_ID* sct_id{nullptr};
+ ATH_CHECK(detStore()->retrieve(sct_id, "SCT_ID"));
+ m_InducedChargeModel = std::make_unique<SCT_InducedChargeModel>(sct_id->wafer_hash_max());
+ }
+
// Surface drift time calculation Stuff
m_tHalfwayDrift = m_tSurfaceDrift * 0.5;
m_distHalfInterStrip = m_distInterStrip * 0.5;
@@ -304,7 +314,7 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
const SiHit& phit,
const ISiSurfaceChargesInserter& inserter,
float p_eventTime,
- unsigned short p_eventId, CLHEP::HepRandomEngine * rndmEngine) const {
+ unsigned short p_eventId, CLHEP::HepRandomEngine* rndmEngine) const {
const SCT_ModuleSideDesign* design{dynamic_cast<const SCT_ModuleSideDesign*>(&(element->design()))};
if (design==nullptr) {
ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::process can not get " << design);
@@ -353,6 +363,24 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
float yhit{xPhi};
float zhit{xDep};
+ SCT_InducedChargeModel::SCT_InducedChargeModelData* data{nullptr};
+ if (m_doInducedChargeModel) { // Setting magnetic field for the ICM.
+ SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, Gaudi::Hive::currentContext()};
+ const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+ float vdepl{m_vdepl};
+ float vbias{m_vbias};
+ if (m_useSiCondDB) {
+ vdepl = m_siConditionsTool->depletionVoltage(hashId);
+ vbias = m_siConditionsTool->biasVoltage(hashId);
+ }
+ data = m_InducedChargeModel->setWaferData(vdepl,
+ vbias,
+ element,
+ fieldCondObj,
+ m_siConditionsTool,
+ rndmEngine);
+ }
+
if (m_doDistortions) {
if (element->isBarrel()) {// Only apply disortions to barrel modules
Amg::Vector2D BOW;
@@ -420,8 +448,12 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
if (t_drift > 0.0) {
const float x1{xhit + StepX * dstep};
float y1{yhit + StepY * dstep};
- y1 += tanLorentz * zReadout; // !< Taking into account the magnetic field
- const float sigma{diffusionSigma(zReadout, element)};
+
+ float sigma{0.};
+ if (not m_doInducedChargeModel) {
+ sigma = diffusionSigma(zReadout, element);
+ y1 += tanLorentz * zReadout; // !< Taking into account the magnetic field
+ } // These are treated in Induced Charge Model.
for (int i{0}; i < m_numberOfCharges; ++i) {
const float rx{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
@@ -491,21 +523,59 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
} // m_doRamo==true
} // chargeIsTrapped()
} // m_doTrapping==true
-
+
if (not m_doRamo) {
- const SiLocalPosition position{element->hitLocalToLocal(xd, yd)};
- if (design->inActiveArea(position)) {
- const float sdist{static_cast<float>(design->scaledDistanceToNearestDiode(position))}; // !< dist on the surface from the hit point to the nearest strip (diode)
- const float t_surf{surfaceDriftTime(2.0 * sdist)}; // !< Surface drift time
- const float totaltime{(m_tfix > -998.) ? m_tfix.value() : t_drift + timeOfFlight + t_surf}; // !< Total drift time
- inserter(SiSurfaceCharge(position, SiCharge(q1, totaltime, hitproc, trklink)));
- } else {
- ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) << "Local position (phi, eta, depth): ("
- << position.xPhi() << ", " << position.xEta() << ", " << position.xDepth()
- << ") of the element is out of active area, charge = " << q1);
+ if (m_doInducedChargeModel) { // Induced Charge Model
+ // Charges storages for 50 ns. 0.5 ns steps.
+ double Q_m2[SCT_InducedChargeModel::NTransportSteps]={0};
+ double Q_m1[SCT_InducedChargeModel::NTransportSteps]={0};
+ double Q_00[SCT_InducedChargeModel::NTransportSteps]={0};
+ double Q_p1[SCT_InducedChargeModel::NTransportSteps]={0};
+ double Q_p2[SCT_InducedChargeModel::NTransportSteps]={0};
+
+ const double mm2cm = 0.1; // For mm -> cm conversion
+ // Unit for y and z : mm -> cm in SCT_InducedChargeModel
+ m_InducedChargeModel->holeTransport(*data,
+ y0*mm2cm, z0*mm2cm,
+ Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+ hashId, m_siPropertiesTool);
+ m_InducedChargeModel->electronTransport(*data,
+ y0*mm2cm, z0*mm2cm,
+ Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+ hashId, m_siPropertiesTool);
+
+ for (int it{0}; it<SCT_InducedChargeModel::NTransportSteps; it++) {
+ if (Q_00[it] == 0.0) continue;
+ double ICM_time{(it+0.5)*0.5 + timeOfFlight};
+ double Q_new[SCT_InducedChargeModel::NStrips]{
+ Q_m2[it], Q_m1[it], Q_00[it], Q_p1[it], Q_p2[it]
+ };
+ for (int strip{SCT_InducedChargeModel::StartStrip}; strip<=SCT_InducedChargeModel::EndStrip; strip++) {
+ double ystrip{y1 + strip * stripPitch};
+ SiLocalPosition position{element->hitLocalToLocal(x1, ystrip)};
+ if (design->inActiveArea(position)) {
+ inserter(SiSurfaceCharge(position,
+ SiCharge(q1 * Q_new[strip+SCT_InducedChargeModel::Offset],
+ ICM_time, hitproc, trklink)));
+ }
+ }
+ }
+ } else { // not m_doInducedChargeModel
+ const SiLocalPosition position{element->hitLocalToLocal(xd, yd)};
+ if (design->inActiveArea(position)) {
+ const float sdist{static_cast<float>(design->scaledDistanceToNearestDiode(position))};
+ // !< dist on the surface from the hit point to the nearest strip (diode)
+ const float t_surf{surfaceDriftTime(2.0 * sdist)}; // !< Surface drift time
+ const float totaltime{(m_tfix > -998.) ? m_tfix.value() : t_drift + timeOfFlight + t_surf}; // !< Total drift time
+ inserter(SiSurfaceCharge(position, SiCharge(q1, totaltime, hitproc, trklink)));
+ } else {
+ ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) << "Local position (phi, eta, depth): ("
+ << position.xPhi() << ", " << position.xEta() << ", " << position.xDepth()
+ << ") of the element is out of active area, charge = " << q1);
+ }
}
- } // end of loop on charges
- }
+ }
+ } // end of loop on charges
}
}
return;
diff --git a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
index a6a89fabcdc0d2883aa31a48fe56b084ac48b865..9510fc80b29c4f616a061e8aa04787fe7b86a19c 100644
--- a/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
+++ b/InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
@@ -1,7 +1,7 @@
// -*- C++ -*-
/*
- Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
/**
@@ -36,12 +36,16 @@
#include "AthenaBaseComps/AthAlgTool.h"
#include "SCT_Digitization/ISCT_SurfaceChargesGenerator.h"
+#include "SCT_InducedChargeModel.h"
+
#include "Identifier/IdentifierHash.h"
#include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
#include "InDetCondTools/ISiLorentzAngleTool.h"
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
#include "SCT_ConditionsTools/ISCT_RadDamageSummaryTool.h"
#include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
#include "SiPropertiesTool/ISiPropertiesTool.h"
+#include "StoreGate/ReadCondHandle.h"
// Gaudi
#include "GaudiKernel/ITHistSvc.h" // for ITHistSvc
@@ -123,6 +127,7 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
BooleanProperty m_doHistoTrap{this, "doHistoTrap", false, "Histogram the charge trapping effect"};
BooleanProperty m_doRamo{this, "doRamo", false, "Ramo Potential for charge trapping effect"};
BooleanProperty m_isOverlay{this, "isOverlay", false, "flag for overlay"};
+ BooleanProperty m_doInducedChargeModel{this, "doInducedChargeModel", false, "Flag for Induced Charge Model"};
//ToolHandles
ToolHandle<ISCT_ModuleDistortionsTool> m_distortionsTool{this, "SCTDistortionsTool", "SCT_DistortionsTool", "Tool to retrieve SCT distortions"};
@@ -133,6 +138,9 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
ServiceHandle<ITHistSvc> m_thistSvc{this, "THistSvc", "THistSvc"};
+ SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCacheCondObjInputKey{
+ this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
+
float m_tHalfwayDrift{0.}; //!< Surface drift time
float m_distInterStrip{1.0}; //!< Inter strip distance normalized to 1
float m_distHalfInterStrip{0.}; //!< Half way distance inter strip
@@ -162,6 +170,9 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
TProfile* m_h_velocity_trap{nullptr};
TProfile* m_h_mobility_trap{nullptr};
TH1F* m_h_trap_pos{nullptr};
+
+ // For Induced Charge Module, M.Togawa
+ std::unique_ptr<SCT_InducedChargeModel> m_InducedChargeModel;
};
#endif // SCT_SURFACECHARGESGENERATOR_H