23 #include "CLHEP/Geometry/Point3D.h"
24 #include "CLHEP/Random/RandomEngine.h"
25 #include "CLHEP/Random/RandGaussZiggurat.h"
27 #include "GaudiKernel/ConcurrencyFlags.h"
41 const std::string&
name,
69 ATH_MSG_FATAL(
"Filling histograms not supported in MT jobs.");
70 return StatusCode::FAILURE;
75 m_h = std::make_unique<Hists>();
82 const SCT_ID* sct_id{
nullptr};
94 ATH_MSG_INFO(
"\tsurface drift still not on, wrong params");
99 ATH_MSG_FATAL(
"\tCannot set both FixedTime and SubtractTime options!");
100 ATH_MSG_INFO(
"\tMake sure the two flags are not set simultaneously in jo");
101 return StatusCode::FAILURE;
106 return StatusCode::SUCCESS;
178 return StatusCode::SUCCESS;
187 return StatusCode::SUCCESS;
194 if (element==
nullptr) {
195 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::process element is nullptr");
199 if (design==
nullptr) {
200 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::process can not get " << design);
203 const double thickness{design->thickness()};
206 if ((zhit < 0.0) or (zhit > thickness)) {
211 float depletionVoltage{0.};
212 float biasVoltage{0.};
221 const float denominator{
static_cast<float>(depletionVoltage + biasVoltage - (2.0 * zhit * depletionVoltage / thickness))};
223 if (biasVoltage >= depletionVoltage) {
225 ATH_MSG_ERROR(
"driftTime: negative argument X for log(X) " << zhit);
235 t_drift *= thickness * thickness / (2.0 *
m_siPropertiesTool->getSiProperties(hashId, ctx).holeDriftMobility() * depletionVoltage);
243 if (element==
nullptr) {
244 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::diffusionSigma element is nullptr");
251 const float sigma{
static_cast<float>(std::sqrt(2. *
m_siPropertiesTool->getSiProperties(hashId, ctx).holeDiffusionConstant() *
t))};
263 const float sensorThickness{
static_cast<float>(element->
thickness())};
264 return driftTime(sensorThickness, element, ctx);
276 const float sensorThickness{
static_cast<float>(element->
thickness())};
291 float t_surfaceDrift{0.};
299 return t_surfaceDrift;
316 CLHEP::HepRandomEngine * rndmEngine,
317 const EventContext& ctx) {
330 unsigned short p_eventId,
331 CLHEP::HepRandomEngine* rndmEngine,
332 const EventContext& ctx)
const {
335 if (initialDesign==
nullptr) {
336 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::process can not get " << initialDesign);
346 if(motherDesign!=
nullptr){
348 design = motherDesign;
351 ATH_MSG_DEBUG(
"No Mother Design - Using Design from DetElement directly!");
352 design = initialDesign;
355 const double thickness{design->
thickness()};
363 float timeOfFlight{p_eventTime +
hitTime(phit)};
381 const float cEta{
static_cast<float>(endPos[
SiHit::xEta]) - xEta};
382 const float cPhi{
static_cast<float>(endPos[
SiHit::xPhi]) - xPhi};
383 const float cDep{
static_cast<float>(endPos[
SiHit::xDep]) - xDep};
387 const float largeStep{std::sqrt(cEta*cEta + cPhi*cPhi + cDep*cDep)};
391 const float q1{
static_cast<float>(
e1 *
m_siPropertiesTool->getSiProperties(hashId, ctx).electronHolePairsPerEnergy())};
427 const float stepX{cX / numberOfSteps};
428 const float stepY{cY / numberOfSteps};
429 const float stepZ{cZ / numberOfSteps};
446 for (
int istep{0}; istep < numberOfSteps; ++istep) {
448 float z1{zhit + stepZ * dstep};
452 float zReadout{
static_cast<float>(0.5 * thickness - design->
readoutSide() * z1)};
453 const double spess{zReadout};
457 h.m_h_depD->Fill(z1);
458 h.m_h_spess->Fill(spess);
461 float t_drift{
driftTime(zReadout, element, ctx)};
462 if (t_drift>-2.0000002 and t_drift<-1.9999998) {
463 ATH_MSG_DEBUG(
"Checking for rounding errors in compression");
464 if ((std::abs(z1) - 0.5 * thickness) < 0.000010) {
465 ATH_MSG_DEBUG(
"Rounding error found attempting to correct it. z1 = " << std::fixed << std::setprecision(8) << z1);
467 z1 = 0.0000005 - 0.5 * thickness;
470 z1 = 0.5 * thickness - 0.0000005;
473 zReadout = 0.5 * thickness - design->
readoutSide() * z1;
474 t_drift =
driftTime(zReadout, element, ctx);
475 if (t_drift>-2.0000002 and t_drift<-1.9999998) {
478 ATH_MSG_DEBUG(
"Correction Successful! z1 = " << std::fixed << std::setprecision(8) << z1 <<
", zReadout = " << zReadout <<
", t_drift = " << t_drift);
481 ATH_MSG_DEBUG(
"No rounding error found. Making no correction.");
485 const float x1{xhit + stepX * dstep};
486 float y1{yhit + stepY * dstep};
491 y1 += tanLorentz * zReadout;
495 const float rx{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
497 const float ry{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
498 const float yd{
y1 +
sigma * ry};
501 const double stripPitch{0.080};
502 double dstrip{
y1 / stripPitch};
504 dstrip = dstrip - std::trunc(dstrip);
506 dstrip = dstrip - std::trunc(dstrip) + 1;
510 double y0{dstrip * stripPitch};
511 double z0{thickness - zReadout};
517 h.m_h_zhit->Fill(zhit);
519 double trap_pos{-999999.}, drift_time{-999999.};
524 double Q_m2{0.}, Q_m1{0.}, Q_00{0.}, Q_p1{0.}, Q_p2{0.};
526 dstrip =
y1 / stripPitch;
532 dstrip -=
static_cast<double>(
static_cast<int>(dstrip));
534 dstrip -=
static_cast<double>(
static_cast<int>(dstrip)) + 1;
538 double yfin{dstrip * stripPitch};
539 double zfin{thickness - trap_pos};
541 m_radDamageTool->holeTransport(y0,
z0, yfin, zfin, Q_m2, Q_m1, Q_00, Q_p1, Q_p2);
542 for (
int strip{-2}; strip<=2; strip++) {
543 const double ystrip{yd + strip * stripPitch};
547 if (strip == -2)
charge = Q_m2;
548 else if (strip == -1)
charge = Q_m1;
549 else if (strip == 0)
charge = Q_00;
550 else if (strip == 1)
charge = Q_p1;
551 else if (strip == 2)
charge = Q_p2;
552 const double time{drift_time};
573 const double mm2cm = 0.1;
577 Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
582 Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
587 if (Q_00[
it] == 0.0)
continue;
588 double ICM_time{(
it+0.5)*0.5 + timeOfFlight};
590 Q_m2[
it], Q_m1[
it], Q_00[
it], Q_p1[
it], Q_p2[
it]
593 double ystrip{
y1 + strip * stripPitch};
598 ICM_time, hitproc, trklink)));
608 const float totaltime{(
m_tfix > -998.) ?
m_tfix.value() : t_drift + timeOfFlight + t_surf};
611 ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) <<
"Local position (phi, eta, depth): ("
612 << position.xPhi() <<
", " << position.xEta() <<
", " << position.xDepth()
613 <<
") of the element is out of active area, charge = " << q1);
625 double& drift_time)
const {
626 if (element==
nullptr) {
627 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::chargeIsTrapped element is nullptr");
630 bool isTrapped{
false};
633 const double electric_field{condData.getElectricField()};
637 const double mobChar{condData.getHoleDriftMobility()};
638 h.m_h_efieldz->Fill(spess, electric_field);
639 h.m_h_efield->Fill(electric_field);
640 h.m_h_mob_Char->Fill(electric_field, mobChar);
641 h.m_h_vel->Fill(electric_field, electric_field * mobChar);
643 const double t_electrode{condData.getTimeToElectrode()};
644 drift_time = condData.getTrappingTime();
645 const double z_trap{condData.getTrappingPositionZ()};
646 trap_pos = spess - z_trap;
649 h.m_h_drift_time->Fill(drift_time);
650 h.m_h_t_electrode->Fill(t_electrode);
651 h.m_h_drift_electrode->Fill(drift_time, t_electrode);
652 h.m_h_ztrap_tot->Fill(z_trap);
656 if (drift_time < t_electrode) {
658 ATH_MSG_INFO(
"drift_time: " << drift_time <<
" t_electrode: " << t_electrode <<
" spess " << spess);
662 h.m_h_ztrap->Fill(z_trap);
663 h.m_h_trap_drift_t->Fill(drift_time);
664 h.m_h_drift1->Fill(spess, drift_time / t_electrode);
665 h.m_h_gen->Fill(spess, drift_time);
666 h.m_h_gen1->Fill(spess, z_trap);
667 h.m_h_gen2->Fill(spess, z_trap / drift_time * t_electrode);
668 h.m_h_velocity_trap->Fill(electric_field, z_trap / drift_time);
669 h.m_h_mobility_trap->Fill(electric_field, z_trap / drift_time / electric_field);
670 h.m_h_trap_pos->Fill(trap_pos);
675 const double z_trap{condData.getTrappingPositionZ()};
677 h.m_h_no_ztrap->Fill(z_trap);
678 h.m_h_notrap_drift_t->Fill(drift_time);