332 {
333 const SCT_ModuleSideDesign* design;
334 const SCT_ModuleSideDesign* initialDesign{
dynamic_cast<const SCT_ModuleSideDesign*
>(&(element->
design()))};
335 if (initialDesign==nullptr) {
336 ATH_MSG_ERROR(
"StripSurfaceChargesGenerator::process can not get " << initialDesign);
337 return;
338 }
339
340
341
342
343
344 const SCT_ModuleSideDesign* motherDesign = initialDesign->
getMother();
345
346 if(motherDesign!=nullptr){
348 design = motherDesign;
349 }
350 else {
351 ATH_MSG_DEBUG(
"No Mother Design - Using Design from DetElement directly!");
352 design = initialDesign;
353 }
354
355 const double thickness{design->
thickness()};
358
359
360
361
362
363 float timeOfFlight{p_eventTime +
hitTime(phit)};
364
365
366 timeOfFlight -= (element->
center().
mag()) / CLHEP::c_light;
367
368
369
372 }
373
374
379
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};
384
385
386
387 const float largeStep{std::sqrt(cEta*cEta + cPhi*cPhi + cDep*cDep)};
390 const float e1{
static_cast<float>(phit.
energyLoss() / steps)};
391 const float q1{
static_cast<float>(
e1 *
m_siPropertiesTool->getSiProperties(hashId, ctx).electronHolePairsPerEnergy())};
392
393
394
395
396
397
398
399
400
401 float xhit{xDep};
402 float yhit{xPhi};
403 float zhit{xEta};
404 float cX{cDep};
405 float cY{cPhi};
406 float cZ{cEta};
407
408 InducedChargeModel::SCT_InducedChargeModelData*
data{
nullptr};
411 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
414 return;
415 }
421 }
423 vbias,
424 element,
425 fieldCondObj,
427 rndmEngine,
428 ctx);
429 }
430
431 const float stepX{cX / numberOfSteps};
432 const float stepY{cY / numberOfSteps};
433 const float stepZ{cZ / numberOfSteps};
434
435
436
441
442
443
444
446 }
447 }
448
449 float dstep{-0.5};
450 for (int istep{0}; istep < numberOfSteps; ++istep) {
451 dstep += 1.0;
452 float z1{zhit + stepZ * dstep};
453
454
455
456 float zReadout{
static_cast<float>(0.5 * thickness - design->
readoutSide() * z1)};
457 const double spess{zReadout};
458
461 h.m_h_depD->Fill(z1);
462 h.m_h_spess->Fill(spess);
463 }
464
465 float t_drift{
driftTime(zReadout, element, ctx)};
466 if (t_drift>-2.0000002 and t_drift<-1.9999998) {
467 ATH_MSG_DEBUG(
"Checking for rounding errors in compression");
468 if ((std::abs(z1) - 0.5 * thickness) < 0.000010) {
469 ATH_MSG_DEBUG(
"Rounding error found attempting to correct it. z1 = " << std::fixed << std::setprecision(8) << z1);
470 if (z1 < 0.0) {
471 z1 = 0.0000005 - 0.5 * thickness;
472
473 } else {
474 z1 = 0.5 * thickness - 0.0000005;
475
476 }
477 zReadout = 0.5 * thickness - design->
readoutSide() * z1;
478 t_drift =
driftTime(zReadout, element, ctx);
479 if (t_drift>-2.0000002 and t_drift<-1.9999998) {
481 } else {
482 ATH_MSG_DEBUG(
"Correction Successful! z1 = " << std::fixed << std::setprecision(8) << z1 <<
", zReadout = " << zReadout <<
", t_drift = " << t_drift);
483 }
484 } else {
485 ATH_MSG_DEBUG(
"No rounding error found. Making no correction.");
486 }
487 }
488 if (t_drift > 0.0) {
489 const float x1{xhit + stepX * dstep};
490 float y1{yhit + stepY * dstep};
491
495 y1 += tanLorentz * zReadout;
496 }
497
499 const float rx{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
501 const float ry{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
502 const float yd{
y1 +
sigma * ry};
503
504
505 const double stripPitch{0.080};
506 double dstrip{
y1 / stripPitch};
507 if (dstrip > 0.) {
508 dstrip = dstrip - std::trunc(dstrip);
509 } else {
510 dstrip = dstrip - std::trunc(dstrip) + 1;
511 }
512
513
514 double y0{dstrip * stripPitch};
515 double z0{thickness - zReadout};
516
517
521 h.m_h_zhit->Fill(zhit);
522 }
523 double trap_pos{-999999.}, drift_time{-999999.};
526 break;
527 } else {
528 double Q_m2{0.}, Q_m1{0.}, Q_00{0.}, Q_p1{0.}, Q_p2{0.};
529
530 dstrip =
y1 / stripPitch;
531
532
533
534
535 if (dstrip > 0.) {
536 dstrip -= static_cast<double>(static_cast<int>(dstrip));
537 } else {
538 dstrip -= static_cast<double>(static_cast<int>(dstrip)) + 1;
539 }
540
541
542 double yfin{dstrip * stripPitch};
543 double zfin{thickness - trap_pos};
544
545 m_radDamageTool->holeTransport(y0, z0, yfin, zfin, Q_m2, Q_m1, Q_00, Q_p1, Q_p2, ctx);
547 const double ystrip{yd +
strip * stripPitch};
556 const double time{drift_time};
558 inserter(SiSurfaceCharge(position, SiCharge(q1*
charge, time, hitproc, trklink)));
559 continue;
560 }
561 }
562 }
564 }
565 }
566 }
567
570
576
577 const double mm2cm = 0.1;
578
580 y0*mm2cm, z0*mm2cm,
581 Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
583 ctx);
585 y0*mm2cm, z0*mm2cm,
586 Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
588 ctx);
589
591 if (Q_00[it] == 0.0) continue;
592 double ICM_time{(
it+0.5)*0.5 + timeOfFlight};
594 Q_m2[
it], Q_m1[
it], Q_00[
it], Q_p1[
it], Q_p2[
it]
595 };
597 double ystrip{
y1 +
strip * stripPitch};
600 inserter(SiSurfaceCharge(position,
602 ICM_time, hitproc, trklink)));
603 }
604 }
605 }
606 } else {
610
612 const float totaltime{(
m_tfix > -998.) ?
m_tfix.value() : t_drift + timeOfFlight + t_surf};
613 inserter(SiSurfaceCharge(position, SiCharge(q1, totaltime, hitproc, trklink)));
614 } else {
615 ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) <<
"Local position (phi, eta, depth): ("
617 << ") of the element is out of active area, charge = " << q1);
618 }
619 }
620 }
621 }
622 }
623 }
624 }
float hitTime(const AFP_SIDSimHit &hit)
Scalar mag() const
mag method
#define ATH_MSG_WARNING(x)
double charge(const T &p)
bool isValid() const
Validity check.
static HepMcParticleLink getRedirectedLink(const HepMcParticleLink &particleLink, uint32_t eventIndex, const EventContext &ctx)
Return a HepMcParticleLink pointing at the same particle, but in a different GenEvent.
float surfaceDriftTime(float ysurf) const
Calculate of the surface drift time.
IntegerProperty m_numberOfCharges
bool chargeIsTrapped(double spess, const InDetDD::SiDetectorElement *element, double &trap_pos, double &drift_time, const EventContext &ctx) const
FloatProperty m_smallStepLength
int readoutSide() const
ReadoutSide.
virtual double scaledDistanceToNearestDiode(const SiLocalPosition &chargePos) const =0
give distance to the nearest diode in units of pitch, from 0.0 to 0.5, this method should be fast as ...
const SCT_ModuleSideDesign * getMother() const
virtual bool inActiveArea(const SiLocalPosition &chargePos, bool checkBondGap=true) const =0
check if the position is in active area
virtual const Amg::Vector3D & center() const override final
Center in global coordinates.
Amg::Vector2D hitLocalToLocal(double xEta, double xPhi) const
Simulation/Hit local frame to reconstruction local frame.
double energyLoss() const
HepGeom::Point3D< double > localStartPosition() const
const HepMcParticleLink & particleLink() const
HepGeom::Point3D< double > localEndPosition() const
time(flags, cells_name, *args, **kw)
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
double e1(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 1st sampling