diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx index ad5a14bf9908754a1df41ce181383734bfd5ab8e..a7c40e14f060119f77ed854ad2d5c464da30d204 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx @@ -144,8 +144,15 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr& phit, const PixelHistoConverter& avgChargeMap_e = fluenceData->getAvgChargeMap3D_e(); const PixelHistoConverter& avgChargeMap_h = fluenceData->getAvgChargeMap3D_h(); - for (unsigned int istep = 0; istep < trfHitRecord.size(); istep++) { - std::pair < double, double > iHitRecord = trfHitRecord[istep]; + //Parameters which will be smeared by a random value for each charge propagated, + //recomputed for every trfHitRecord but initialized only once + std::vector DtElectron (ncharges, 0.); + std::vector DtHole (ncharges, 0.); + std::vector rdifElectron (ncharges, 0.); + std::vector rdifHole (ncharges, 0.); + + for (size_t istep = 0; istep < trfHitRecord.size(); istep++) { + std::pair< double,double> const & iHitRecord = trfHitRecord[istep]; double eta_i = eta_0; double phi_i = phi_0; @@ -182,9 +189,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr& phit, int nPixX = int(x_new / pixel_size_x); int nPixY = int(y_new / pixel_size_y); ATH_MSG_DEBUG(" -- nPixX = " << nPixX << " nPixY = " << nPixY); - // In case the charge moves into a neighboring pixel - int extraNPixX = nPixX; - int extraNPixY = nPixY; //position relative to the bottom left corner of the pixel double x_pix = x_new - pixel_size_x * (nPixX); @@ -206,117 +210,290 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr& phit, const double mobilityElectron = getMobility(efield, false); const double mobilityHole = getMobility(efield, true); - + auto driftTimeElectron = getDriftTime(false, ncharges); + auto driftTimeHole = getDriftTime(true, ncharges); //Need to determine how many elementary charges this charge chunk represents. double chunk_size = energy_per_step * eleholePairEnergy; //number of electrons/holes - //set minimum limit to prevent dividing into smaller subcharges than one fundamental charge if (chunk_size < 1) chunk_size = 1; const double kappa = 1. / std::sqrt(chunk_size); + + const double timeToElectrodeElectron = timeMap_e.getContent(timeMap_e.getBinX(1e3*y_pix), timeMap_e.getBinY(1e3*x_pix)); + const double timeToElectrodeHole = timeMap_h.getContent(timeMap_h.getBinX(1e3*y_pix), timeMap_h.getBinY(1e3*x_pix)); + + /* Diffusion via the Einstein relation + D = mu * kB * T / q + D = (mu / mm^2/MV*ns) * (T/273 K) * 0.024 microns^2 / ns */ + const auto prefactor_e = mobilityElectron*0.024*m_temperature / 273.; + const auto prefactor_h = mobilityHole*0.024*m_temperature / 273.; + //Apply diffusion. rdif is the max. diffusion + for(size_t i = 0; i < ncharges; ++i) { + DtElectron[i] = prefactor_e * std::min(driftTimeElectron[i], timeToElectrodeElectron); + DtHole[i] = prefactor_h * std::min(driftTimeHole[i], timeToElectrodeHole); + rdifElectron[i] = 1e-3*std::sqrt(DtElectron[i]); //in mm + rdifHole[i] = 1e-3*std::sqrt(DtHole[i]); //in mm + } - //Loop over charge-carrier pairs + const float average_chargeElectron = avgChargeMap_e.getContent(avgChargeMap_e.getBinY(1e3*y_pix), avgChargeMap_e.getBinX(1e3*x_pix)); + const float average_chargeHole = avgChargeMap_h.getContent(avgChargeMap_h.getBinY(1e3*y_pix), avgChargeMap_h.getBinX(1e3*x_pix)); + + //We're sticking to the "old" convention here where there was a loop over -1/0/1 for both directions + //(x/y) to the neighbouring pixels. We call them (p)lus, (z)ero, (m)inus in either i- or j- direction + //hence giving us plus-in-i or pi or zero-in-j or zj etc. + //i is in direction of x, whereas j is into the y direction, this code is ugly, but given that we don't have very few branching conditions + //it's pretty fast + + //This was: ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3 - i * pixel_size_x) + const std::size_t ramo_init_bin_y_pi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 2)); + const std::size_t ramo_init_bin_y_zi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3)); + const std::size_t ramo_init_bin_y_mi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 4)); + + //This was: ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y - j * pixel_size_y) + const std::size_t ramo_init_bin_x_pj = ramoPotentialMap.getBinX(1000*(y_pix - 0.5*pixel_size_y)); + const std::size_t ramo_init_bin_x_zj = ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y)); + const std::size_t ramo_init_bin_x_mj = ramoPotentialMap.getBinX(1000*(y_pix + 1.5*pixel_size_y)); + + //For some reason the x- and y-indices for some values are swapped, this is extremely confusing + float ramoInit_pipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_pi); + float ramoInit_pizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_pi); + float ramoInit_pimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_pi); + float ramoInit_zipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_zi); + float ramoInit_zizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_zi); + float ramoInit_zimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_zi); + float ramoInit_mipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_mi); + float ramoInit_mizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_mi); + float ramoInit_mimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_mi); + + auto mc_particle_link = HepMcParticleLink(phit->trackNumber(), phit.eventId(), evColl, idxFlag, ctx); + const auto hit_time = hitTime(phit); + + //Loop over charge-carrier pairs, we're looping over electrons and holes at the same time for (int j = 0; j < ncharges; j++) { - // Loop over everything twice: once for electrons and once for holes - for (int eholes = 0; eholes < 2; eholes++) { - const bool isHole = (eholes == 1); // Set a condition to keep track of electron/hole-specific functions - - // Reset extraPixel coordinates each time through loop - extraNPixX = nPixX; - extraNPixY = nPixY; - - double timeToElectrode(0.0); //[ns] - if (!isHole) { - timeToElectrode = timeMap_e.getContent(timeMap_e.getBinX(1e3*y_pix),timeMap_e.getBinY(1e3*x_pix)); - } - else { - timeToElectrode = timeMap_h.getContent(timeMap_h.getBinX(1e3*y_pix),timeMap_h.getBinY(1e3*x_pix)); - } - - const double driftTime = getDriftTime(isHole); - - //Apply drift due to diffusion - const double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine); - - //Apply diffusion. rdif is teh max. diffusion - const double Dt = (isHole ? mobilityHole : mobilityElectron) * (0.024) * std::min(driftTime, timeToElectrode) * m_temperature / 273.; - const double rdif = 1e-3*std::sqrt(Dt); //in mm - double xposDiff = x_pix + rdif * phiRand; - const double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine); - double yposDiff = y_pix + rdif * etaRand; + // In case the charge moves into a neighboring pixel + int extraNPixXElectron = nPixX; + int extraNPixYElectron = nPixY; + int extraNPixXHole = nPixX; + int extraNPixYHole = nPixY; + + //Apply drift due to diffusion + std::array randomNumbers; + CLHEP::RandGaussZiggurat::shootArray(rndmEngine, 4, randomNumbers.data()); + + double xposDiffElectron = x_pix + rdifElectron[j] * randomNumbers[0]; + double yposDiffElectron = y_pix + rdifElectron[j] * randomNumbers[1]; + double xposDiffHole = x_pix + rdifHole[j] * randomNumbers[2]; + double yposDiffHole = y_pix + rdifHole[j] * randomNumbers[3]; // Account for drifting into another pixel - while (xposDiff > pixel_size_x) { - extraNPixX = extraNPixX + 1; // increments or decrements pixel count in x - xposDiff = xposDiff - pixel_size_x; // moves xpos coordinate 1 pixel over in x + while (xposDiffElectron > pixel_size_x) { + extraNPixXElectron++; // increments or decrements pixel count in x + xposDiffElectron -= pixel_size_x; // moves xpos coordinate 1 pixel over in x } - while (xposDiff < 0) { - extraNPixX = extraNPixX - 1; - xposDiff = xposDiff + pixel_size_x; + while (xposDiffElectron < 0) { + extraNPixXElectron--; + xposDiffElectron += pixel_size_x; } - while (yposDiff > pixel_size_y) { - extraNPixY = extraNPixY + 1; // increments or decrements pixel count in y - yposDiff = yposDiff - pixel_size_y; // moves xpos coordinate 1 pixel over in y + while (yposDiffElectron > pixel_size_y) { + extraNPixYElectron++; // increments or decrements pixel count in y + yposDiffElectron -= pixel_size_y; // moves xpos coordinate 1 pixel over in y } - while (yposDiff < 0) { - extraNPixY = extraNPixY - 1; - yposDiff = yposDiff + pixel_size_y; + while (yposDiffElectron < 0) { + extraNPixYElectron--; + yposDiffElectron += pixel_size_y; } - float average_charge = isHole ? avgChargeMap_h.getContent(avgChargeMap_h.getBinY(1e3*y_pix),avgChargeMap_h.getBinX(1e3*x_pix)) : - avgChargeMap_e.getContent(avgChargeMap_e.getBinY(1e3*y_pix),avgChargeMap_e.getBinX(1e3*x_pix)); - - double xposFinal(0); // Trapping position Y - double yposFinal(0); // Trapping position X - if (!isHole) { - xposFinal = yPositionMap_e.getContent(yPositionMap_e.getBinX(1e3*yposDiff),yPositionMap_e.getBinY(1e3*xposDiff),yPositionMap_e.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm] - yposFinal = xPositionMap_e.getContent(xPositionMap_e.getBinX(1e3*yposDiff),xPositionMap_e.getBinY(1e3*xposDiff),xPositionMap_e.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm] - } - else { - xposFinal = yPositionMap_h.getContent(yPositionMap_h.getBinX(1e3*yposDiff),yPositionMap_h.getBinY(1e3*xposDiff),yPositionMap_h.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm] - yposFinal = xPositionMap_h.getContent(xPositionMap_h.getBinX(1e3*yposDiff),xPositionMap_h.getBinY(1e3*xposDiff),xPositionMap_h.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm] + //And drifting for for holes + while (xposDiffHole > pixel_size_x) { + extraNPixXHole++; + xposDiffHole -= pixel_size_x; + } + while (xposDiffHole < 0) { + extraNPixXHole--; + xposDiffHole += pixel_size_x; + } + while (yposDiffHole > pixel_size_y) { + extraNPixYHole++; + yposDiffHole -= pixel_size_y; + } + while (yposDiffHole < 0) { + extraNPixYHole--; + yposDiffHole += pixel_size_y; } + auto xposFinalElectron = yPositionMap_e.getContent( yPositionMap_e.getBinX(1e3*yposDiffElectron), + yPositionMap_e.getBinY(1e3*xposDiffElectron), + yPositionMap_e.getBinZ(std::min(driftTimeElectron[j],timeToElectrodeElectron)) )*1e-3; //[mm] + auto yposFinalElectron = xPositionMap_e.getContent( xPositionMap_e.getBinX(1e3*yposDiffElectron), + xPositionMap_e.getBinY(1e3*xposDiffElectron), + xPositionMap_e.getBinZ(std::min(driftTimeElectron[j],timeToElectrodeElectron)) )*1e-3; //[mm] + + auto xposFinalHole = yPositionMap_h.getContent( yPositionMap_h.getBinX(1e3*yposDiffHole), + yPositionMap_h.getBinY(1e3*xposDiffHole), + yPositionMap_h.getBinZ(std::min(driftTimeHole[j],timeToElectrodeHole)) )*1e-3; //[mm] + auto yposFinalHole = xPositionMap_h.getContent( xPositionMap_h.getBinX(1e3*yposDiffHole), + xPositionMap_h.getBinY(1e3*xposDiffHole), + xPositionMap_h.getBinZ(std::min(driftTimeHole[j],timeToElectrodeHole)) )*1e-3; //[mm] // -- Calculate signal in current pixel and in the neighboring ones - // -- loop in the x-coordinate - for (int i = -1; i <= 1; i++) { - double xNeighbor = i * pixel_size_x; - // -- loop in the y-coordinate - const std::size_t ramo_init_bin_y = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3 - xNeighbor)); - const std::size_t ramo_final_bin_y = ramoPotentialMap.getBinY(1000*(xposFinal + pixel_size_x * 3 - xNeighbor)); - - for (int j = -1; j <= 1; j++) { - double yNeighbor = j * pixel_size_y; - - //Ramo map over 500umx350um pixel area - //Ramo init different if charge diffused into neighboring pixel -> change primary pixel!! - float ramoInit = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y); - float ramoFinal = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(1000*(yposFinal + 0.5*pixel_size_y - yNeighbor)), ramo_final_bin_y); - - // Record deposit - double eHitRamo = (isHole ? -1. : 1.) * energy_per_step * (ramoFinal - ramoInit); - - if (doChunkCorrection) { - eHitRamo = energy_per_step * average_charge + kappa * (eHitRamo - energy_per_step * average_charge); - } - - double induced_charge = eHitRamo * eleholePairEnergy; - - // -- pixel coordinates --> module coordinates - 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); - - const SiSurfaceCharge scharge(chargePos, - SiCharge(induced_charge, - pHitTime, SiCharge::track, particleLink)); - const SiCellId& diode = Module.cellIdOfPosition(scharge.position()); - if (diode.isValid()) { - const SiCharge& charge = scharge.charge(); - chargedDiodes.add(diode, charge); - } - } - } - } + //This was: ramoPotentialMap.getBinY(1000*(xposFinal + pixel_size_x * 3 - i * pixel_size_x) + const std::size_t ramo_final_bin_y_pi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 2)); + const std::size_t ramo_final_bin_y_pi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 2)); + const std::size_t ramo_final_bin_y_zi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 3)); + const std::size_t ramo_final_bin_y_zi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 3)); + const std::size_t ramo_final_bin_y_mi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 4)); + const std::size_t ramo_final_bin_y_mi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 4)); + //This was: ramoPotentialMap.getBinX(1000*(yposFinal + 0.5*pixel_size_y - j * pixel_size_y) + const std::size_t ramo_final_bin_x_pj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron - 0.5*pixel_size_y)); + const std::size_t ramo_final_bin_x_pj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole - 0.5*pixel_size_y)); + const std::size_t ramo_final_bin_x_zj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron + 0.5*pixel_size_y)); + const std::size_t ramo_final_bin_x_zj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole + 0.5*pixel_size_y)); + const std::size_t ramo_final_bin_x_mj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron + 1.5*pixel_size_y)); + const std::size_t ramo_final_bin_x_mj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole + 1.5*pixel_size_y)); + + float ramoFinal_pipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_pi_electrons); + float ramoFinal_pizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_pi_electrons); + float ramoFinal_pimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_pi_electrons); + float ramoFinal_zipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_zi_electrons); + float ramoFinal_zizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_zi_electrons); + float ramoFinal_zimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_zi_electrons); + float ramoFinal_mipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_mi_electrons); + float ramoFinal_mizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_mi_electrons); + float ramoFinal_mimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_mi_electrons); + + float ramoFinal_pipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_pi_holes); + float ramoFinal_pizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_pi_holes); + float ramoFinal_pimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_pi_holes); + float ramoFinal_zipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_zi_holes); + float ramoFinal_zizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_zi_holes); + float ramoFinal_zimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_zi_holes); + float ramoFinal_mipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_mi_holes); + float ramoFinal_mizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_mi_holes); + float ramoFinal_mimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_mi_holes); + + // Record deposit + double eHitRamo_pipj_electrons = energy_per_step * (ramoFinal_pipj_electrons - ramoInit_pipj); + double eHitRamo_pipj_holes = -1. * energy_per_step * (ramoFinal_pipj_holes - ramoInit_pipj); + double eHitRamo_pizj_electrons = energy_per_step * (ramoFinal_pizj_electrons - ramoInit_pizj); + double eHitRamo_pizj_holes = -1. * energy_per_step * (ramoFinal_pizj_holes - ramoInit_pizj); + double eHitRamo_pimj_electrons = energy_per_step * (ramoFinal_pimj_electrons - ramoInit_pimj); + double eHitRamo_pimj_holes = -1. * energy_per_step * (ramoFinal_pimj_holes - ramoInit_pimj); + double eHitRamo_zipj_electrons = energy_per_step * (ramoFinal_zipj_electrons - ramoInit_zipj); + double eHitRamo_zipj_holes = -1. * energy_per_step * (ramoFinal_zipj_holes - ramoInit_zipj); + double eHitRamo_zizj_electrons = energy_per_step * (ramoFinal_zizj_electrons - ramoInit_zizj); + double eHitRamo_zizj_holes = -1. * energy_per_step * (ramoFinal_zizj_holes - ramoInit_zizj); + double eHitRamo_zimj_electrons = energy_per_step * (ramoFinal_zimj_electrons - ramoInit_zimj); + double eHitRamo_zimj_holes = -1. * energy_per_step * (ramoFinal_zimj_holes - ramoInit_zimj); + double eHitRamo_mipj_electrons = energy_per_step * (ramoFinal_mipj_electrons - ramoInit_mipj); + double eHitRamo_mipj_holes = -1. * energy_per_step * (ramoFinal_mipj_holes - ramoInit_mipj); + double eHitRamo_mizj_electrons = energy_per_step * (ramoFinal_mizj_electrons - ramoInit_mizj); + double eHitRamo_mizj_holes = -1. * energy_per_step * (ramoFinal_mizj_holes - ramoInit_mizj); + double eHitRamo_mimj_electrons = energy_per_step * (ramoFinal_mimj_electrons - ramoInit_mimj); + double eHitRamo_mimj_holes = -1. * energy_per_step * (ramoFinal_mimj_holes - ramoInit_mimj); + + if(doChunkCorrection) { + eHitRamo_pipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pipj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_pipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pipj_holes - energy_per_step * average_chargeHole); + eHitRamo_pizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pizj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_pizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pizj_holes - energy_per_step * average_chargeHole); + eHitRamo_pimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pimj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_pimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pimj_holes - energy_per_step * average_chargeHole); + eHitRamo_zipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zipj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_zipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_zipj_holes - energy_per_step * average_chargeHole); + eHitRamo_zizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zizj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_zizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_zizj_holes - energy_per_step * average_chargeHole); + eHitRamo_zimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zimj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_pimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pimj_holes - energy_per_step * average_chargeHole); + eHitRamo_mipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mipj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_mipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mipj_holes - energy_per_step * average_chargeHole); + eHitRamo_mizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mizj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_mizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mizj_holes - energy_per_step * average_chargeHole); + eHitRamo_mimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mimj_electrons - energy_per_step * average_chargeElectron); + eHitRamo_mimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mimj_holes - energy_per_step * average_chargeHole); + } + + // -- pixel coordinates --> module coordinates + //This was: double x_mod = x_pix + pixel_size_x * extraNPixX - 0.5*module_size_x + i*pixel_size_x; + //This was: double y_mod = y_pix + pixel_size_y * extraNPixY - 0.5*module_size_y + j*pixel_size_y; + double x_mod_pi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x + pixel_size_x; + double y_mod_pj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y + pixel_size_y; + double x_mod_pi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x + pixel_size_x; + double y_mod_pj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y + pixel_size_y; + double x_mod_zi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x; + double y_mod_zj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y; + double x_mod_zi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x; + double y_mod_zj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y; + double x_mod_mi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x - pixel_size_x; + double y_mod_mj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y - pixel_size_y; + double x_mod_mi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x - pixel_size_x; + double y_mod_mj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y - pixel_size_y; + + //This was: const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod) (for whatever reason half the maps x/y are inverted...) + const SiLocalPosition& chargePos_pipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_pi_electrons); + const SiLocalPosition& chargePos_pizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_pi_electrons); + const SiLocalPosition& chargePos_pimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_pi_electrons); + const SiLocalPosition& chargePos_zipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_zi_electrons); + const SiLocalPosition& chargePos_zizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_zi_electrons); + const SiLocalPosition& chargePos_zimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_zi_electrons); + const SiLocalPosition& chargePos_mipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_mi_electrons); + const SiLocalPosition& chargePos_mizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_mi_electrons); + const SiLocalPosition& chargePos_mimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_mi_electrons); + + const SiLocalPosition& chargePos_pipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_pi_holes); + const SiLocalPosition& chargePos_pizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_pi_holes); + const SiLocalPosition& chargePos_pimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_pi_holes); + const SiLocalPosition& chargePos_zipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_zi_holes); + const SiLocalPosition& chargePos_zizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_zi_holes); + const SiLocalPosition& chargePos_zimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_zi_holes); + const SiLocalPosition& chargePos_mipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_mi_holes); + const SiLocalPosition& chargePos_mizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_mi_holes); + const SiLocalPosition& chargePos_mimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_mi_holes); + + //This was: const SiSurfaceCharge scharge(chargePos, SiCharge(eHitRamo*eleholePairEnergy, hit_time, SiCharge::track, particleLink)); + const SiSurfaceCharge scharge_pipj_electrons(chargePos_pipj_electrons, SiCharge(eHitRamo_pipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_pizj_electrons(chargePos_pizj_electrons, SiCharge(eHitRamo_pizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_pimj_electrons(chargePos_pimj_electrons, SiCharge(eHitRamo_pimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zipj_electrons(chargePos_zipj_electrons, SiCharge(eHitRamo_zipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zizj_electrons(chargePos_zizj_electrons, SiCharge(eHitRamo_zizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zimj_electrons(chargePos_zimj_electrons, SiCharge(eHitRamo_zimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mipj_electrons(chargePos_mipj_electrons, SiCharge(eHitRamo_mipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mizj_electrons(chargePos_mizj_electrons, SiCharge(eHitRamo_mizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mimj_electrons(chargePos_mimj_electrons, SiCharge(eHitRamo_mimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + + const SiSurfaceCharge scharge_pipj_holes(chargePos_pipj_holes, SiCharge(eHitRamo_pipj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_pizj_holes(chargePos_pizj_holes, SiCharge(eHitRamo_pizj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_pimj_holes(chargePos_pimj_holes, SiCharge(eHitRamo_pimj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zipj_holes(chargePos_zipj_holes, SiCharge(eHitRamo_zipj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zizj_holes(chargePos_zizj_holes, SiCharge(eHitRamo_zizj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_zimj_holes(chargePos_zimj_holes, SiCharge(eHitRamo_zimj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mipj_holes(chargePos_mipj_holes, SiCharge(eHitRamo_mipj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mizj_holes(chargePos_mizj_holes, SiCharge(eHitRamo_mizj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + const SiSurfaceCharge scharge_mimj_holes(chargePos_mimj_holes, SiCharge(eHitRamo_mimj_holes*eleholePairEnergy, hit_time, SiCharge::track, mc_particle_link)); + + auto addCharge = [&Module, &chargedDiodes](SiSurfaceCharge const & scharge) { + const SiCellId& diode = Module.cellIdOfPosition(scharge.position()); + if (diode.isValid()) { + const SiCharge& charge = scharge.charge(); + chargedDiodes.add(diode, charge); + } + }; + + addCharge(scharge_pipj_electrons); + addCharge(scharge_pizj_electrons); + addCharge(scharge_pimj_electrons); + addCharge(scharge_zipj_electrons); + addCharge(scharge_zizj_electrons); + addCharge(scharge_zimj_electrons); + addCharge(scharge_mipj_electrons); + addCharge(scharge_mizj_electrons); + addCharge(scharge_mimj_electrons); + addCharge(scharge_pipj_holes); + addCharge(scharge_pizj_holes); + addCharge(scharge_pimj_holes); + addCharge(scharge_zipj_holes); + addCharge(scharge_zizj_holes); + addCharge(scharge_zimj_holes); + addCharge(scharge_mipj_holes); + addCharge(scharge_mizj_holes); + addCharge(scharge_mimj_holes); } } } else { @@ -516,17 +693,18 @@ double SensorSim3DTool::getMobility(double electricField, bool isHoleBit) { return mobility; // mm^2/(MV*ns) } -double SensorSim3DTool::getDriftTime(bool isHoleBit) { - double u = CLHEP::RandFlat::shoot(0., 1.); // - double driftTime = 0; - - // need to update to std::logf when we update gcc - this is a known bug in gcc libc - if (isHoleBit) { - driftTime = (-1.) * m_trappingTimeHoles * logf(u); // ns - } else { - driftTime = (-1.) * m_trappingTimeElectrons * logf(u); // ns +std::vector SensorSim3DTool::getDriftTime(bool isHoleBit, size_t n) { + std::vector rand (n, 0.); + std::vector result (n, 0.); + CLHEP::RandFlat::shootArray(n, rand.data(), 0., 1.); + for(size_t i = 0; i < n; i++) { + if (isHoleBit) { + result[i]= (-1.) * m_trappingTimeHoles * logf(rand[i]); // ns + } else { + result[i] = (-1.) * m_trappingTimeElectrons * logf(rand[i]); // ns + } } - return driftTime; + return result; } diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h index 88713a2593ff971704b451ad663fb00e24273fc1..b4d06d18aa5fde912430c25c5b1eccac3a77a346 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h @@ -42,7 +42,7 @@ public: StatusCode printProbMap(const std::string&) const; double getMobility(double electricField, bool isHoleBit); - double getDriftTime(bool isHoleBit); + std::vector getDriftTime(bool isHoleBit, size_t number); private: