diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
index ef5ebb4c0c6ac4b1912e5789db9b1481f0a9a22c..1bbb33b98eef2bba39409a7ad941e487700718a0 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
@@ -287,8 +287,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
chargepos.z());
// -- change origin of coordinates to the left bottom of module
- double x_new = chargepos.x() + module_size_x / 2.;
- double y_new = chargepos.y() + module_size_y / 2.;
+ double x_new = chargepos.x() + 0.5*module_size_x;
+ double y_new = chargepos.y() + 0.5*module_size_y;
// -- change from module frame to pixel frame
@@ -322,11 +322,10 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
double eHit = energy_per_step;
//Need to determine how many elementary charges this charge chunk represents.
double chunk_size = energy_per_step * eleholePairEnergy; //number of electrons/holes
- ATH_MSG_DEBUG("Chunk size: " << energy_per_step << "*" << eleholePairEnergy << " = " << chunk_size);
//set minimum limit to prevent dividing into smaller subcharges than one fundamental charge
if (chunk_size < 1) chunk_size = 1;
- double kappa = 1. / sqrt(chunk_size);
+ double kappa = 1. / std::sqrt(chunk_size);
// Loop over everything twice: once for electrons and once for holes
for (int eholes = 0; eholes < 2; eholes++) {
@@ -346,7 +345,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//Apply diffusion. rdif is teh max. diffusion
double Dt =
getMobility(efield, isHole) * (0.024) * std::min(driftTime, timeToElectrode) * m_temperature / 273.;
- double rdif = sqrt(Dt) / 1000; //in mm
+ double rdif = 1e-3*std::sqrt(Dt); //in mm
double xposDiff = x_pix + rdif * phiRand;
double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
double yposDiff = y_pix + rdif * etaRand;
@@ -370,18 +369,12 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
}
- ATH_MSG_DEBUG(" -- diffused position w.r.t. pixel edge = " << xposDiff << " " << yposDiff);
-
float average_charge = isHole ? m_avgChargeMap_h.getContent(m_avgChargeMap_h.getBinY(1e3*y_pix), m_avgChargeMap_h.getBinX(1e3*x_pix)) :
m_avgChargeMap_e.getContent(m_avgChargeMap_e.getBinY(1e3*y_pix), m_avgChargeMap_e.getBinX(1e3*x_pix));
- ATH_MSG_DEBUG(" -- driftTime, timeToElectrode = " << driftTime << " " << timeToElectrode);
-
double xposFinal = getTrappingPositionY(yposDiff, xposDiff, std::min(driftTime, timeToElectrode), isHole);
double yposFinal = getTrappingPositionX(yposDiff, xposDiff, std::min(driftTime, timeToElectrode), isHole);
- ATH_MSG_DEBUG(" -- trapped position w.r.t. pixel edge = " << xposFinal << " " << yposFinal);
-
// -- Calculate signal in current pixel and in the neighboring ones
// -- loop in the x-coordinate
for (int i = -1; i <= 1; i++) {
@@ -393,13 +386,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
for (int j = -1; j <= 1; j++) {
double yNeighbor = j * pixel_size_y;
- ATH_MSG_DEBUG(
- " -- Ramo init position w.r.t. Ramo map edge = " << x_pix + pixel_size_x * 3 - xNeighbor << " " << y_pix + pixel_size_y * 1 / 2 -
- yNeighbor);
- ATH_MSG_DEBUG(
- " -- Ramo final position w.r.t. Ramo map edge = " << xposFinal + pixel_size_x * 3 - xNeighbor << " " << yposFinal + pixel_size_y * 1 / 2 -
- yNeighbor);
-
//Ramo map over 500umx350um pixel area
//Ramo init different if charge diffused into neighboring pixel -> change primary pixel!!
float ramoInit = m_ramoPotentialMap[index].getContent(m_ramoPotentialMap[index].getBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y);
@@ -408,31 +394,24 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
// Record deposit
double eHitRamo = (1 - 2 * isHole) * eHit * (ramoFinal - ramoInit);
- ATH_MSG_DEBUG(
- "At neighbor pixel " << i << " " << j << " Hit of " << eHitRamo << " including Ramo factor: " << ramoFinal -
- ramoInit);
-
if (m_doChunkCorrection) {
- ATH_MSG_DEBUG("Energy before chunk correction: " << eHitRamo);
eHitRamo = eHit * average_charge + kappa * (eHitRamo - eHit * average_charge);
- ATH_MSG_DEBUG("Energy after chunk correction: " << eHitRamo);
}
double induced_charge = eHitRamo * eleholePairEnergy;
// -- pixel coordinates --> module coordinates
- double x_mod = x_pix + xNeighbor + pixel_size_x * extraNPixX - module_size_x / 2.;
- double y_mod = y_pix + yNeighbor + pixel_size_y * extraNPixY - module_size_y / 2.;
- SiLocalPosition chargePos = Module.hitLocalToLocal(y_mod, x_mod);
+ double x_mod = x_pix + xNeighbor + pixel_size_x * extraNPixX - 0.5*module_size_x;
+ double y_mod = y_pix + yNeighbor + pixel_size_y * extraNPixY - 0.5*module_size_y;
+ const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod);
- SiSurfaceCharge scharge(chargePos, SiCharge(induced_charge, hitTime(
+ const SiSurfaceCharge scharge(chargePos, SiCharge(induced_charge, hitTime(
phit), SiCharge::track, HepMcParticleLink(
phit->trackNumber(), phit.eventId(), evColl, idxFlag, ctx)));
- SiCellId diode = Module.cellIdOfPosition(scharge.position());
- SiCharge charge = scharge.charge();
+ const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
if (diode.isValid()) {
+ const SiCharge& charge = scharge.charge();
chargedDiodes.add(diode, charge);
- ATH_MSG_DEBUG("induced charge: " << induced_charge << " x_mod: " << x_mod << " y_mod: " << y_mod);
}
}
}
@@ -474,8 +453,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
chargepos.z());
// -- change origin of coordinates to the left bottom of module
- double x_new = chargepos.x() + module_size_x / 2.;
- double y_new = chargepos.y() + module_size_y / 2.;
+ double x_new = chargepos.x() + 0.5*module_size_x;
+ double y_new = chargepos.y() + 0.5*module_size_y;
// -- change from module frame to pixel frame
int nPixX = int(x_new / pixel_size_x);
@@ -500,32 +479,32 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
y_neighbor = y_pix_center - j * pixel_size_y;
// -- check if the neighbor falls inside the charge collection prob map window
- if ((fabs(x_neighbor) < pixel_size_x) && (fabs(y_neighbor) < pixel_size_y)) {
+ if ((std::abs(x_neighbor) < pixel_size_x) && (std::abs(y_neighbor) < pixel_size_y)) {
// -- change origin of coordinates to the bottom left of the charge
// collection prob map "window", i.e. shift of 1-pixel twd bottom left
double x_neighbor_map = x_neighbor + pixel_size_x;
double y_neighbor_map = y_neighbor + pixel_size_y;
- int x_bin_cc_map = static_cast<int>(x_neighbor_map / x_bin_size);
- int y_bin_cc_map = static_cast<int>(y_neighbor_map / y_bin_size);
+ int x_bin_cc_map = x_neighbor_map / x_bin_size;
+ int y_bin_cc_map = y_neighbor_map / y_bin_size;
// -- retrieve the charge collection probability from Svc
// -- swap x and y bins to match Map coord convention
- double ccprob_neighbor = getProbMapEntry("FEI4", y_bin_cc_map, x_bin_cc_map);
+ double ccprob_neighbor = getProbMapEntry(SensorType::FEI4, y_bin_cc_map, x_bin_cc_map);
if (ccprob_neighbor == -1.) return StatusCode::FAILURE;
double ed = es_current * eleholePairEnergy * ccprob_neighbor;
// -- pixel coordinates --> module coordinates
- double x_mod = x_neighbor + pixel_size_x / 2 + pixel_size_x * nPixX - module_size_x / 2.;
- double y_mod = y_neighbor + pixel_size_y / 2 + pixel_size_y * nPixY - module_size_y / 2.;
- SiLocalPosition chargePos = Module.hitLocalToLocal(y_mod, x_mod);
+ double x_mod = x_neighbor + 0.5*pixel_size_x + pixel_size_x * nPixX - 0.5*module_size_x;
+ double y_mod = y_neighbor + 0.5*pixel_size_y + pixel_size_y * nPixY - 0.5*module_size_y;
+ const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod);
- SiSurfaceCharge scharge(chargePos, SiCharge(ed, hitTime(phit), SiCharge::track, HepMcParticleLink(
+ const SiSurfaceCharge scharge(chargePos, SiCharge(ed, hitTime(phit), SiCharge::track, HepMcParticleLink(
phit->trackNumber(), phit.eventId(), evColl, idxFlag, ctx)));
- SiCellId diode = Module.cellIdOfPosition(scharge.position());
- SiCharge charge = scharge.charge();
+ const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
if (diode.isValid()) {
+ const SiCharge& charge = scharge.charge();
chargedDiodes.add(diode, charge);
}
}
@@ -595,13 +574,13 @@ StatusCode SensorSim3DTool::printProbMap(const std::string& readout) const {
}
// -- Returns the Charge Collection Probability at a given point (bin_x,bin_y)
-double SensorSim3DTool::getProbMapEntry(const std::string& readout, int binx, int biny) const {
+double SensorSim3DTool::getProbMapEntry(const SensorType& readout, int binx, int biny) const {
std::pair<int, int> doublekey(binx, biny);
double echarge;
- if (readout == "FEI4") {
+ if (readout == SensorType::FEI4) {
std::multimap<std::pair<int, int>, double>::const_iterator iter = m_probMapFEI4.find(doublekey);
echarge = iter->second;
- } else if (readout == "FEI3") {
+ } else if (readout == SensorType::FEI3) {
std::multimap<std::pair<int, int>, double>::const_iterator iter = m_probMapFEI3.find(doublekey);
echarge = iter->second;
} else {
@@ -630,17 +609,17 @@ double SensorSim3DTool::getMobility(double electricField, bool isHoleBit) {
// https://cds.cern.ch/record/684187/files/indet-2001-004.pdf).
if (!isHoleBit) {
- vsat = 15.3 * pow(m_temperature, -0.87); // mm/ns
- ecrit = 1.01E-7 * pow(m_temperature, 1.55); // MV/mm
- beta = 2.57E-2 * pow(m_temperature, 0.66);
+ vsat = 15.3 * std::pow(m_temperature, -0.87); // mm/ns
+ ecrit = 1.01E-7 * std::pow(m_temperature, 1.55); // MV/mm
+ beta = 2.57E-2 * std::pow(m_temperature, 0.66);
}
if (isHoleBit) {
- vsat = 1.62 * pow(m_temperature, -0.52); // mm/ns
- ecrit = 1.24E-7 * pow(m_temperature, 1.68); // MV/mm
- beta = 0.46 * pow(m_temperature, 0.17);
+ vsat = 1.62 * std::pow(m_temperature, -0.52); // mm/ns
+ ecrit = 1.24E-7 * std::pow(m_temperature, 1.68); // MV/mm
+ beta = 0.46 * std::pow(m_temperature, 0.17);
}
- double mobility = (vsat / ecrit) / pow(1 + pow((electricField / ecrit), beta), (1 / beta));
+ double mobility = (vsat / ecrit) / std::pow(1 + std::pow((electricField / ecrit), beta), (1 / beta));
return mobility; // mm^2/(MV*ns)
}
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
index ea37734ebf756fcdfd8c4b0c37f0c7fc8b5a0ba7..91bda4541b9b5570dac0fdb9176e0c0a6aca98c4 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
@@ -41,7 +41,6 @@ public:
// 3D sensor simulation using probability density map (used in RUN-2 (no radiation damage)
StatusCode readProbMap(const std::string&);
StatusCode printProbMap(const std::string&) const;
- double getProbMapEntry(const std::string&, int, int) const;
double getElectricField(double x, double y);
double getMobility(double electricField, bool isHoleBit);
@@ -50,7 +49,15 @@ public:
double getTrappingPositionX(double initX, double initY, double driftTime, bool isHoleBit);
double getTrappingPositionY(double initX, double initY, double driftTime, bool isHoleBit);
private:
+
+ enum class SensorType {
+ FEI4,
+ FEI3
+ };
+
SensorSim3DTool();
+
+ double getProbMapEntry(const SensorType&, int, int) const;
// 3D sensor simulation using probability density map (used in RUN-2 (no radiation damage)
std::multimap<std::pair<int, int>, double> m_probMapFEI4;
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
index d0df06e475e3c309643d528fbfcaa0ee76a61e1f..bd72e4a0130c8b6e6d77a68ed2e2425863a4bab1 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
@@ -245,7 +245,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
double collectionDist = 0.2 * CLHEP::mm;
double smearScale = 1. + 0.35 * smearRand;
double tanLorentz = m_lorentzAngleTool->getTanLorentzAngle(Module.identifyHash());
- double coLorentz = std::sqrt(1.0 + pow(tanLorentz, 2));
+ double coLorentz = std::sqrt(1.0 + (tanLorentz*tanLorentz));
const EBC_EVCOLL evColl = EBC_MAINEVCOLL;
const HepMcParticleLink::PositionFlag idxFlag =
@@ -304,56 +304,52 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
nontrappingProbability = exp(-dist_electrode / collectionDist);
}
- const std::size_t distance_f_e_bin_x = m_doInterpolateEfield ? m_distanceMap_e[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_e(layer).getBinX(dist_electrode);
- const std::size_t distance_f_h_bin_x = m_doInterpolateEfield ? m_distanceMap_h[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_h(layer).getBinX(dist_electrode);
- const std::size_t tanLorentz_e_bin_x = m_doInterpolateEfield ? m_lorentzMap_e[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_e(layer).getBinX(dist_electrode);
- const std::size_t tanLorentz_h_bin_x = m_doInterpolateEfield ? m_lorentzMap_h[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_h(layer).getBinX(dist_electrode);
+ if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
+
+ const PixelHistoConverter& distanceMap_e = m_doInterpolateEfield ? m_distanceMap_e[layer] : moduleData->getDistanceMap_e(layer);
+ const PixelHistoConverter& distanceMap_h = m_doInterpolateEfield ? m_distanceMap_h[layer] : moduleData->getDistanceMap_h(layer);
+ const PixelHistoConverter& lorentzMap_e = m_doInterpolateEfield ? m_lorentzMap_e[layer] : moduleData->getLorentzMap_e(layer);
+ const PixelHistoConverter& lorentzMap_h = m_doInterpolateEfield ? m_lorentzMap_h[layer] : moduleData->getLorentzMap_h(layer);
+ const PixelHistoConverter& ramoPotentialMap = m_doInterpolateEfield ? m_ramoPotentialMap[layer] : moduleData->getRamoPotentialMap(layer);
+
+ const std::size_t distance_f_e_bin_x = distanceMap_e.getBinX(dist_electrode);
+ const std::size_t distance_f_h_bin_x = distanceMap_h.getBinX(dist_electrode);
+ const std::size_t tanLorentz_e_bin_x = lorentzMap_e.getBinX(dist_electrode);
+ const std::size_t tanLorentz_h_bin_x = lorentzMap_h.getBinX(dist_electrode);
- for (int j = 0; j < ncharges; j++) {
- if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
+ for (int j = 0; j < ncharges; j++) {
double u = CLHEP::RandFlat::shoot(0., 1.);
- double drifttime_e = (-1.) * (trappingTimes.first) * TMath::Log(u); //ns
+ const double drifttime_e = (-1.) * (trappingTimes.first) * std::log(u); //ns
u = CLHEP::RandFlat::shoot(0., 1.);
- double drifttime_h = (-1.) * (trappingTimes.second) * TMath::Log(u); //ns
+ const double drifttime_h = (-1.) * (trappingTimes.second) * std::log(u); //ns
//Now, need the z-position at the trap.
- double depth_f_e = 0.0;
- double depth_f_h = 0.0;
- double tanLorentz_e = 0.0;
- double tanLorentz_h = 0.0;
//TODO: the holes map does not currently extend for a drift time long enough that, any hole will reach
//the corresponding electrode. This needs to be rectified by either (a) extrapolating the current map or
//(b) making a new map with a y-axis (drift time) that extends to at least 18 ns so all charge carriers reach
// electrode.
//However, if choose (b), will need to reduce granularity of map.
- if (m_doInterpolateEfield) {
- depth_f_e = m_distanceMap_e[layer].getContent(distance_f_e_bin_x, m_distanceMap_e[layer].getBinY(drifttime_e));
- depth_f_h = m_distanceMap_h[layer].getContent(distance_f_h_bin_x, m_distanceMap_h[layer].getBinY(drifttime_h));
- tanLorentz_e = m_lorentzMap_e[layer].getContent(tanLorentz_e_bin_x, m_lorentzMap_e[layer].getBinY(depth_f_e));
- tanLorentz_h = m_lorentzMap_h[layer].getContent(tanLorentz_h_bin_x, m_lorentzMap_h[layer].getBinY(depth_f_h));
- } else { // use fluence value from conditions data
- depth_f_e = moduleData->getDistanceMap_e(layer).getContent(distance_f_e_bin_x, moduleData->getDistanceMap_e(layer).getBinY(drifttime_e));
- depth_f_h = moduleData->getDistanceMap_h(layer).getContent(distance_f_h_bin_x, moduleData->getDistanceMap_h(layer).getBinY(drifttime_h));
- tanLorentz_e = moduleData->getLorentzMap_e(layer).getContent(tanLorentz_e_bin_x, moduleData->getLorentzMap_e(layer).getBinY(depth_f_e));
- tanLorentz_h = moduleData->getLorentzMap_h(layer).getContent(tanLorentz_h_bin_x, moduleData->getLorentzMap_h(layer).getBinY(depth_f_h));
- }
- double dz_e = fabs(dist_electrode - depth_f_e);
- double dz_h = fabs(depth_f_h - dist_electrode);
- double coLorentz_e = std::sqrt(1.0 + std::pow(tanLorentz_e, 2));
+ const double depth_f_e = distanceMap_e.getContent(distance_f_e_bin_x, distanceMap_e.getBinY(drifttime_e));
+ const double depth_f_h = distanceMap_h.getContent(distance_f_h_bin_x, distanceMap_h.getBinY(drifttime_h));
+ const double tanLorentz_e = lorentzMap_e.getContent(tanLorentz_e_bin_x, lorentzMap_e.getBinY(depth_f_e));
+ const double tanLorentz_h = lorentzMap_h.getContent(tanLorentz_h_bin_x, lorentzMap_h.getBinY(depth_f_h));
+ const double dz_e = std::abs(dist_electrode - depth_f_e);
+ const double dz_h = std::abs(depth_f_h - dist_electrode);
+ const double coLorentz_e = std::sqrt(1.0 + (tanLorentz_e*tanLorentz_e));
//Apply drift due to Lorentz force and diffusion
double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
//Apply diffusion. rdif is teh max. diffusion
- double rdif_e = this->m_diffusionConstant * sqrt(fabs(dist_electrode - depth_f_e) * coLorentz_e / 0.3);
- double phi_f_e = phi_i + dz_e * tanLorentz_e + rdif_e * phiRand;
+ const double rdif_e = this->m_diffusionConstant * std::sqrt(std::abs(dist_electrode - depth_f_e) * coLorentz_e / 0.3);
+ const double phi_f_e = phi_i + dz_e * tanLorentz_e + rdif_e * phiRand;
double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
double eta_f_e = eta_i + rdif_e * etaRand;
phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
- double coLorentz_h = std::sqrt(1.0 + std::pow(tanLorentz_h, 2));
- double rdif_h = this->m_diffusionConstant * sqrt(fabs(dist_electrode - depth_f_h) * coLorentz_h / 0.3);
- double phi_f_h = phi_i + dz_h * tanLorentz_h + rdif_h * phiRand;
+ const double coLorentz_h = std::sqrt(1.0 + (tanLorentz_h*tanLorentz_h));
+ const double rdif_h = this->m_diffusionConstant * std::sqrt(std::abs(dist_electrode - depth_f_h) * coLorentz_h / 0.3);
+ const double phi_f_h = phi_i + dz_h * tanLorentz_h + rdif_h * phiRand;
etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
double eta_f_h = eta_i + rdif_h * etaRand;
@@ -362,15 +358,15 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
// Slim Edge for IBL planar sensors:
if (p_design.getReadoutTechnology() == InDetDD::PixelModuleDesign::FEI4) {
- ATH_CHECK(applyIBLSlimEdges(energy_per_step, eta_f_e));
- ATH_CHECK(applyIBLSlimEdges(energy_per_step, eta_f_h));
+ applyIBLSlimEdges(energy_per_step, eta_f_e);
+ applyIBLSlimEdges(energy_per_step, eta_f_h);
}
- const std::size_t ramo_f_e_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_e);
- const std::size_t ramo_f_h_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_h);
+ const std::size_t ramo_f_e_bin_z = ramoPotentialMap.getBinZ(1e3*depth_f_e);
+ const std::size_t ramo_f_h_bin_z = ramoPotentialMap.getBinZ(1e3*depth_f_h);
- const bool isFirstZ_e = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isFirstZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).isFirstZ(1e3*depth_f_e);
- const bool isOverflowZ_h = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isOverflowZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).isOverflowZ(1e3*depth_f_h);
+ const bool isFirstZ_e = ramoPotentialMap.isFirstZ(1e3*depth_f_e);
+ const bool isOverflowZ_h = ramoPotentialMap.isOverflowZ(1e3*depth_f_h);
//Loop over nearest neighbours in x and y
//We assume that the lateral diffusion is minimal
@@ -394,8 +390,8 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//What is the displacement of the nn pixel from the primary pixel.
//This is to index the correct entry in the Ramo weighting potential map
- double dPhi_nn_centre = centreOfPixel_nn.xPhi() - centreOfPixel_i.xPhi(); //in mm
- double dEta_nn_centre = centreOfPixel_nn.xEta() - centreOfPixel_i.xEta(); //in mm
+ const double dPhi_nn_centre = centreOfPixel_nn.xPhi() - centreOfPixel_i.xPhi(); //in mm
+ const double dEta_nn_centre = centreOfPixel_nn.xEta() - centreOfPixel_i.xEta(); //in mm
//This all has to be done relative to the (0,0) position since the
//Ramo weighting potential is only mapped out for 1/8th of a pixel. Much of this logic is reflecting the
@@ -403,38 +399,28 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//carrier across the boundaries.
//Find the displacment of the charge carriers from the centre of the pixel in +ve quadrant
- double pixelEta_f_e = eta_f_e - centreOfPixel_i.xEta();
- double pixelPhi_f_e = phi_f_e - centreOfPixel_i.xPhi();
+ const double pixelEta_f_e = eta_f_e - centreOfPixel_i.xEta();
+ const double pixelPhi_f_e = phi_f_e - centreOfPixel_i.xPhi();
- double pixelEta_f_h = eta_f_h - centreOfPixel_i.xEta();
- double pixelPhi_f_h = phi_f_h - centreOfPixel_i.xPhi();
+ const double pixelEta_f_h = eta_f_h - centreOfPixel_i.xEta();
+ const double pixelPhi_f_h = phi_f_h - centreOfPixel_i.xPhi();
//Final position of charge carriers wrt nn centre
- double dEta_f_e = pixelEta_f_e - dEta_nn_centre;
- double dPhi_f_e = pixelPhi_f_e - dPhi_nn_centre;
- dEta_f_e *= scale_f;
- double dEta_f_h = pixelEta_f_h - dEta_nn_centre;
- double dPhi_f_h = pixelPhi_f_h - dPhi_nn_centre;
- dEta_f_h *= scale_f;
+ const double dEta_f_e = std::abs(pixelEta_f_e - dEta_nn_centre)*scale_f;
+ const double dPhi_f_e = std::abs(pixelPhi_f_e - dPhi_nn_centre);
+ const double dEta_f_h = 1e3*std::abs(pixelEta_f_h - dEta_nn_centre)*scale_f;
+ const double dPhi_f_h = 1e3*std::abs(pixelPhi_f_h - dPhi_nn_centre);
//Boundary check on maps
double ramo_f_e = 0.0;
double ramo_f_h = 0.0;
if (!isFirstZ_e) {
- if (m_doInterpolateEfield) {
- ramo_f_e = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_e)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z);
- } else {
- ramo_f_e = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_e)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z);
- }
+ ramo_f_e = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(1e3*dPhi_f_e), ramoPotentialMap.getBinY(1e3*dEta_f_e), ramo_f_e_bin_z);
}
if (!isOverflowZ_h) {
- if (m_doInterpolateEfield) {
- ramo_f_h = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_h)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z);
- } else {
- ramo_f_h = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_h)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z);
- }
+ ramo_f_h = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(dPhi_f_h), ramoPotentialMap.getBinY(dEta_f_h), ramo_f_h_bin_z);
}
//Account for the imperfect binning that would cause charge to be double-counted
if (isOverflowZ_h) {
@@ -442,9 +428,9 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
}
if (isFirstZ_e) {
- if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) {
+ if (dEta_f_e >= 0.5*Module.etaPitch() || dPhi_f_e >= 0.5*Module.phiPitch()) {
ramo_f_e = 0.0;
- } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) {
+ } else if (dEta_f_e < 0.5*Module.etaPitch() && dPhi_f_e < 0.5*Module.phiPitch()) {
ramo_f_e = 1.0;
}
}
@@ -452,27 +438,31 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//Given final position of charge carrier, find induced charge. The difference in Ramo weighting potential
// gives the fraction of charge induced.
//The energy_per_step is transformed into charge with the eleholePair per Energy
- double induced_charge = (ramo_f_e - ramo_f_h) * energy_per_step * eleholePairEnergy;
+ const double potentialDiff = ramo_f_e - ramo_f_h;
+ // this variable ^ can be used to apply some cut to skip the loop
+ const double induced_charge = potentialDiff * energy_per_step * eleholePairEnergy;
//Collect charge in centre of each pixel, since location within pixel doesn't matter for record
- SiLocalPosition chargePos = Module.hitLocalToLocal(centreOfPixel_nn.xEta(), centreOfPixel_nn.xPhi());
+ const SiLocalPosition& chargePos = Module.hitLocalToLocal(centreOfPixel_nn.xEta(), centreOfPixel_nn.xPhi());
//The following lines are adapted from SiDigitization's Inserter class
- SiSurfaceCharge scharge(
+ const SiSurfaceCharge scharge(
chargePos,
SiCharge(induced_charge, hitTime(phit), SiCharge::track, particleLink)
);
- SiCellId diode = Module.cellIdOfPosition(scharge.position());
- SiCharge charge = scharge.charge();
+ const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
if (diode.isValid()) {
+ const SiCharge& charge = scharge.charge();
chargedDiodes.add(diode, charge);
} //IF
} //For q
} //for p
- } else { //If no radDamage, run original
+ }
+ } else { //If no radDamage, run original
+ for (int j = 0; j < ncharges; j++) {
// diffusion sigma
- double rdif = this->m_diffusionConstant * sqrt(dist_electrode * coLorentz / 0.3);
+ double rdif = this->m_diffusionConstant * std::sqrt(dist_electrode * coLorentz / 0.3);
// position at the surface
double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
@@ -485,11 +475,11 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
// Slim Edge for IBL planar sensors:
if (!(Module.isDBM()) && p_design.getReadoutTechnology() == InDetDD::PixelModuleDesign::FEI4) {
- ATH_CHECK(applyIBLSlimEdges(energy_per_step, eta_drifted));
+ applyIBLSlimEdges(energy_per_step, eta_drifted);
}
// Get the charge position in Reconstruction local coordinates.
- SiLocalPosition chargePos = Module.hitLocalToLocal(eta_drifted, phi_drifted);
+ const SiLocalPosition& chargePos = Module.hitLocalToLocal(eta_drifted, phi_drifted);
// The parametrization of the sensor efficiency (if needed)
double ed = 0;
@@ -500,13 +490,12 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
}
//The following lines are adapted from SiDigitization's Inserter class
- SiSurfaceCharge scharge(chargePos, SiCharge(ed, hitTime(phit), SiCharge::track, particleLink));
+ const SiSurfaceCharge scharge(chargePos, SiCharge(ed, hitTime(phit), SiCharge::track, particleLink));
- SiCellId diode = Module.cellIdOfPosition(scharge.position());
-
- SiCharge charge = scharge.charge();
+ const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
if (diode.isValid()) {
+ const SiCharge& charge = scharge.charge();
chargedDiodes.add(diode, charge);
}
} //else: no radDamage, run original
@@ -515,11 +504,11 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
return StatusCode::SUCCESS;
}
-StatusCode SensorSimPlanarTool::applyIBLSlimEdges(double& energy_per_step, double& eta_drifted) {
- if (fabs(eta_drifted) > 20.440) {
+void SensorSimPlanarTool::applyIBLSlimEdges(double& energy_per_step, double& eta_drifted) const{
+ if (std::abs(eta_drifted) > 20.440) {
energy_per_step = 0.0;
}
- if (fabs(eta_drifted) < 20.440 && fabs(eta_drifted) > 20.200) {
+ if (std::abs(eta_drifted) < 20.440 && std::abs(eta_drifted) > 20.200) {
if (eta_drifted > 0) {
energy_per_step = energy_per_step * (68.13 - eta_drifted * 3.333);
eta_drifted = eta_drifted - 0.250;
@@ -528,7 +517,7 @@ StatusCode SensorSimPlanarTool::applyIBLSlimEdges(double& energy_per_step, doubl
eta_drifted = eta_drifted + 0.250;
}
}
- if (fabs(eta_drifted) < 20.200 && fabs(eta_drifted) > 20.100) {
+ if (std::abs(eta_drifted) < 20.200 && std::abs(eta_drifted) > 20.100) {
if (eta_drifted > 0) {
energy_per_step = energy_per_step * (41.2 - eta_drifted * 2.0);
eta_drifted = eta_drifted - 0.250;
@@ -537,5 +526,4 @@ StatusCode SensorSimPlanarTool::applyIBLSlimEdges(double& energy_per_step, doubl
eta_drifted = eta_drifted + 0.250;
}
}
- return StatusCode::SUCCESS;
}
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
index 8dd2757762e66cd44bb101716e65329d98c76e0f..3c4b69e527156218d6abec8a57f62108f0be94b0 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
@@ -36,10 +36,11 @@ public:
CLHEP::HepRandomEngine* rndmEngine,
const EventContext &ctx) override;
- //Apply slim edge inefficiencies for IBL sensors
- StatusCode applyIBLSlimEdges(double& energyPerStep, double& eta_drifted);
private:
SensorSimPlanarTool();
+
+ //Apply slim edge inefficiencies for IBL sensors
+ void applyIBLSlimEdges(double& energyPerStep, double& eta_drifted) const;
// Map for radiation damage simulation
std::vector<PixelHistoConverter> m_ramoPotentialMap;