268 return StatusCode::SUCCESS;
273 if (!Module.
isDBM()) {
274 return StatusCode::SUCCESS;
283 double eta_0 = initialConditions[0];
284 double phi_0 = initialConditions[1];
285 double depth_0 = initialConditions[2];
286 double dEta = initialConditions[3];
287 double dPhi = initialConditions[4];
288 double dDepth = initialConditions[5];
289 double ncharges = initialConditions[6];
290 double iTotalLength = initialConditions[7];
292 const double oneOverNchargesTimes1e6 = 1./(1.E+6 * ncharges);
299 int etaCells = p_design.
columns();
300 int phiCells = p_design.
rows();
302 double eleholePairEnergy = 0;
303 double smearRand = 0;
305 double diffusionConstant;
306 if (Module.
isDBM()) {
307 eleholePairEnergy = 1. / (13. *
CLHEP::eV);
308 diffusionConstant = .00265;
309 smearRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
312 diffusionConstant = .007;
316 double smearScale = 1. + 0.35 * smearRand;
317 double tanLorentz(0);
321 coLorentz = std::sqrt(1.0 + (tanLorentz*tanLorentz));
329 const double pHitTime =
hitTime(phit);
331 const double halfEtaPitch = 0.5*Module.
etaPitch();
332 const double halfPhiPitch = 0.5*Module.
phiPitch();
338 std::pair<double, double> trappingTimes;
352 std::map<unsigned, std::pair<SiLocalPosition, double>> cachedChargeMap;
353 std::map<unsigned, SiCellId> diodeCellMap;
354 for (
const auto& iHitRecord : trfHitRecord) {
356 double eta_i = eta_0;
357 double phi_i = phi_0;
358 double depth_i = depth_0;
360 eta_i += 1.0 * iHitRecord.first / iTotalLength *
dEta;
361 phi_i += 1.0 * iHitRecord.first / iTotalLength *
dPhi;
362 depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
368 if (!pixel_i.
isValid())
continue;
372 double dist_electrode = 0.5 * sensorThickness - Module.
design().
readoutSide() * depth_i;
373 if (dist_electrode < 0) {
381 const int nnLoop_pixelEtaMin =
std::max(-1, pixel_i.
etaIndex() + 1 - etaCells);
383 const int nnLoop_pixelPhiMin =
std::max(-1, pixel_i.
phiIndex() + 1 - phiCells);
385 std::array<double, 3> sensorScales{};
387 const std::size_t distance_f_e_bin_x = distanceMap_e.
getBinX(dist_electrode);
388 const std::size_t distance_f_h_bin_x = distanceMap_h.
getBinX(dist_electrode);
389 const std::size_t tanLorentz_e_bin_x = lorentzMap_e.
getBinX(dist_electrode);
390 const std::size_t tanLorentz_h_bin_x = lorentzMap_h.
getBinX(dist_electrode);
392 const std::size_t sizePhi = std::abs(nnLoop_pixelPhiMax - nnLoop_pixelPhiMin) + 1;
394 const auto pixel_eta = pixel_i.
etaIndex();
395 const auto pixel_phi = pixel_i.
phiIndex();
398 std::array<std::pair<double,double>,9 > centrePixelNNEtaPhi;
399 for (
int p = nnLoop_pixelEtaMin;
p <= nnLoop_pixelEtaMax;
p++) {
400 const std::size_t ieta =
p - nnLoop_pixelEtaMin;
404 if (std::abs(columnWidth - 0.6) < 1
e-9) {
406 }
else if (std::abs(columnWidth - 0.45) < 1
e-9) {
408 }
else if (std::abs(columnWidth - 0.5) < 1
e-9) {
411 sensorScales[ieta] = scale_f;
413 for (
int q = nnLoop_pixelPhiMin;
q <= nnLoop_pixelPhiMax;
q++) {
416 const std::size_t iphi =
q - nnLoop_pixelPhiMin;
417 const std::size_t
index = iphi + ieta*sizePhi;
418 centrePixelNNEtaPhi.at(
index).first = centreOfPixel_nn.
xEta();
419 centrePixelNNEtaPhi[
index].second = centreOfPixel_nn.
xPhi();
424 for (
int j = 0; j < ncharges; j++) {
426 double energy_per_step = oneOverNchargesTimes1e6 * iHitRecord.second;
427 double u = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0);
429 const double drifttime_e = (-1.) * (trappingTimes.first) * logf(
u);
430 u = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0);
431 const double drifttime_h = (-1.) * (trappingTimes.second) * logf(
u);
439 const double depth_f_e = distanceMap_e.
getContent(distance_f_e_bin_x, distanceMap_e.
getBinY(drifttime_e));
440 const double depth_f_h = distanceMap_h.
getContent(distance_f_h_bin_x, distanceMap_h.
getBinY(drifttime_h));
441 const double tanLorentz_e = lorentzMap_e.
getContent(tanLorentz_e_bin_x, lorentzMap_e.
getBinY(depth_f_e));
442 const double tanLorentz_h = lorentzMap_h.
getContent(tanLorentz_h_bin_x, lorentzMap_h.
getBinY(depth_f_h));
443 const double dz_e = std::abs(dist_electrode - depth_f_e);
444 const double dz_h = std::abs(depth_f_h - dist_electrode);
445 const double coLorentz_e = std::sqrt(1.0 + (tanLorentz_e*tanLorentz_e));
448 double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
451 const double rdif_e = diffusionConstant * std::sqrt(dz_e * coLorentz_e / 0.3);
452 const double phi_f_e = phi_i + dz_e * tanLorentz_e + rdif_e * phiRand;
453 double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
454 double eta_f_e = eta_i + rdif_e * etaRand;
456 phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
457 const double coLorentz_h = std::sqrt(1.0 + (tanLorentz_h*tanLorentz_h));
458 const double rdif_h = diffusionConstant * std::sqrt(dz_h * coLorentz_h / 0.3);
459 const double phi_f_h = phi_i + dz_h * tanLorentz_h + rdif_h * phiRand;
460 etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
461 double eta_f_h = eta_i + rdif_h * etaRand;
470 const std::size_t ramo_f_e_bin_z = ramoPotentialMap.
getBinZ(1
e3*depth_f_e);
471 const std::size_t ramo_f_h_bin_z = ramoPotentialMap.
getBinZ(1
e3*depth_f_h);
473 const bool isFirstZ_e = ramoPotentialMap.
isFirstZ(1
e3*depth_f_e);
474 const bool isOverflowZ_h = ramoPotentialMap.
isOverflowZ(1
e3*depth_f_h);
476 const double pixelEta_f_e = eta_f_e - centreOfPixel_i.
xEta();
477 const double pixelPhi_f_e = phi_f_e - centreOfPixel_i.
xPhi();
479 const double pixelEta_f_h = eta_f_h - centreOfPixel_i.
xEta();
480 const double pixelPhi_f_h = phi_f_h - centreOfPixel_i.
xPhi();
484 for (
int p = nnLoop_pixelEtaMin;
p <= nnLoop_pixelEtaMax;
p++) {
485 const std::size_t ieta =
p - nnLoop_pixelEtaMin;
487 for (
int q = nnLoop_pixelPhiMin;
q <= nnLoop_pixelPhiMax;
q++) {
491 const std::size_t iphi =
q - nnLoop_pixelPhiMin;
492 const std::size_t
index = iphi + ieta*sizePhi;
495 const std::pair<double,double>& centrePixelNN = centrePixelNNEtaPhi.at(
index);
496 const double dPhi_nn_centre = centrePixelNN.second - centreOfPixel_i.
xPhi();
497 const double dEta_nn_centre = centrePixelNN.first - centreOfPixel_i.
xEta();
505 const double dEta_f_e = std::abs(pixelEta_f_e - dEta_nn_centre)*sensorScales[ieta];
506 const double dPhi_f_e = std::abs(pixelPhi_f_e - dPhi_nn_centre);
507 const double dEta_f_h = 1
e3*std::abs(pixelEta_f_h - dEta_nn_centre)*sensorScales[ieta];
508 const double dPhi_f_h = 1
e3*std::abs(pixelPhi_f_h - dPhi_nn_centre);
511 double ramo_f_e = 0.0;
512 double ramo_f_h = 0.0;
515 if (dEta_f_e >= halfEtaPitch || dPhi_f_e >= halfPhiPitch) {
528 ramo_f_h = ramoPotentialMap.
getContent(ramoPotentialMap.
getBinX(dPhi_f_h), ramoPotentialMap.
getBinY(dEta_f_h), ramo_f_h_bin_z);
534 const double potentialDiff = ramo_f_e - ramo_f_h;
536 const double induced_charge = potentialDiff * energy_per_step * eleholePairEnergy;
538 unsigned key = (
static_cast<unsigned>(pixel_eta-
p) << 16) |
static_cast<unsigned>(pixel_phi-
q);
539 auto cacheIterator = cachedChargeMap.find(
key);
540 if(cacheIterator == cachedChargeMap.end()) {
541 cachedChargeMap.insert(std::make_pair(
key, std::make_pair(Module.
hitLocalToLocal(centrePixelNN.first, centrePixelNN.second), induced_charge)));
543 cacheIterator->second.second += induced_charge;
550 std::for_each(cachedChargeMap.begin(), cachedChargeMap.end(), [&diodeCellMap, &Module, &chargedDiodes, &pHitTime, &particleLink](
auto& pos_charge_pair){
551 auto& key = pos_charge_pair.first;
552 auto& chargePos = pos_charge_pair.second.first;
553 auto& charge_value = pos_charge_pair.second.second;
555 const SiSurfaceCharge scharge(chargePos, SiCharge(charge_value, pHitTime, SiCharge::track, particleLink));
556 auto diodeIterator = diodeCellMap.find(key);
557 if(diodeIterator == diodeCellMap.end()) diodeIterator = diodeCellMap.insert(std::make_pair(key, Module.cellIdOfPosition(scharge.position()))).first;
558 const SiCellId& thisDiode = diodeIterator->second;
559 if (thisDiode.isValid()) {
560 const SiCharge& charge = scharge.charge();
561 chargedDiodes.add(thisDiode, charge);
575 if (layerToRead > 3) {
583 for (
const auto & iHitRecord : trfHitRecord) {
584 double eta_i = eta_0;
585 double phi_i = phi_0;
586 double depth_i = depth_0;
588 eta_i += 1.0 * iHitRecord.first / iTotalLength *
dEta;
589 phi_i += 1.0 * iHitRecord.first / iTotalLength *
dPhi;
590 depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
598 tanLorentz = lorentzCorrectionHist.
getContent(depthZ);
599 coLorentz = std::sqrt(1.0 + (tanLorentz*tanLorentz));
602 const double chargeCorrection = chargeCorrectionHist.
getContent(depthZ);
605 double nontrappingProbability = 1.0;
606 if (Module.
isDBM()) {
607 nontrappingProbability =
exp(-dist_electrode / collectionDist);
610 for (
int j = 0; j < ncharges; j++) {
612 double energy_per_step = oneOverNchargesTimes1e6 * iHitRecord.second;
614 double rdif = diffusionConstant * std::sqrt(dist_electrode * coLorentz / 0.3);
617 double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
618 double phi_drifted = phi_i + dist_electrode * tanLorentz + rdif * phiRand;
619 double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
620 double eta_drifted = eta_i + rdif * etaRand;
632 if (Module.
isDBM()) {
633 ed = energy_per_step * eleholePairEnergy * nontrappingProbability * smearScale;
635 ed = energy_per_step * eleholePairEnergy;
638 ed *= chargeCorrection;
652 for (
const auto & iHitRecord : trfHitRecord) {
653 double eta_i = eta_0;
654 double phi_i = phi_0;
655 double depth_i = depth_0;
657 eta_i += 1.0 * iHitRecord.first / iTotalLength *
dEta;
658 phi_i += 1.0 * iHitRecord.first / iTotalLength *
dPhi;
659 depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
664 double dist_electrode = 0.5 * sensorThickness - Module.
design().
readoutSide() * depth_i;
665 if (dist_electrode < 0) {
670 double nontrappingProbability = 1.0;
671 if (Module.
isDBM()) {
672 nontrappingProbability =
exp(-dist_electrode / collectionDist);
675 for (
int j = 0; j < ncharges; j++) {
677 double energy_per_step = oneOverNchargesTimes1e6 * iHitRecord.second;
679 double rdif = diffusionConstant * std::sqrt(dist_electrode * coLorentz / 0.3);
682 double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
683 double phi_drifted = phi_i + dist_electrode * tanLorentz + rdif * phiRand;
684 double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
685 double eta_drifted = eta_i + rdif * etaRand;
697 if (Module.
isDBM()) {
698 ed = energy_per_step * eleholePairEnergy * nontrappingProbability * smearScale;
700 ed = energy_per_step * eleholePairEnergy;
716 return StatusCode::SUCCESS;