376{
377
380 rngWrapper->
setSeed( rngName, ctx );
381 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
382
383 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
m_pixelDetEleCollKey, ctx);
384 const InDetDD::SiDetectorElementCollection* elements(*pixelDetEleHandle);
385 if (not pixelDetEleHandle.isValid() or elements==nullptr) {
387 return StatusCode::FAILURE;
388 }
389
391
394
395 SG::ReadCondHandle<PixelChargeCalibCondData> calibDataHandle(
m_chargeDataKey, ctx);
396 const PixelChargeCalibCondData *calibData = *calibDataHandle;
398 std::vector<int> truthIdList;
399 std::vector<Identifier> detEl;
400
402
404
405 truthIdList.clear();
406 detEl.clear();
407
408 while (i != e) {
409
410 const TimedHitPtr<SiHit>&
hit(*i++);
411
412
414 const int barrelEC =
hit->getBarrelEndcap();
415 const int layerDisk =
hit->getLayerDisk();
417 const int etaModule =
hit->getEtaModule();
418
419 const Identifier moduleID =
m_pixel_ID->wafer_id(barrelEC, layerDisk, phiModule, etaModule);
420 const IdentifierHash waferHash =
m_pixel_ID->wafer_hash(moduleID);
421 const InDetDD::SiDetectorElement* hitSiDetElement = elements->getDetectorElement(waferHash);
422 if (!hitSiDetElement) {
ATH_MSG_ERROR(
" could not get detector element ");
continue;}
423
424 if (!(hitSiDetElement->
isPixel())) {
continue;}
425
426
427
428 std::vector<HepMcParticleLink> hit_vector;
429
430 const int truthID = (currentLink.
barcode() !=0 && currentLink.
id() == 0) ? 3 : currentLink.
id();
431
432 const Identifier hitId = hitSiDetElement->
identify();
433
434
435
436 bool isRep = false;
437
438 for (int j : truthIdList) {
439 for (auto & k : detEl) {
440 if ((truthID > 0) && (truthID == j) && (hitId == k)) {isRep = true; break;}
441 }
442 if (isRep) break;
443 }
444
445 if (isRep) continue;
446
447 truthIdList.push_back(truthID);
448 detEl.push_back(hitId);
449
450 HepGeom::Point3D<double> localStartPosition =
hit->localStartPosition();
451 HepGeom::Point3D<double> localEndPosition =
hit->localEndPosition();
452
455
456 int isEndcap = (barrelEC==0) ? 0:1;
457
460
461
462
463
464 Diffuse(localStartPosition, localEndPosition, shiftX, shiftY);
465
466 const Amg::Vector3D localDirection(localEndPosition.x()-localStartPosition.x(), localEndPosition.y()-localStartPosition.y(), localEndPosition.z()-localStartPosition.z());
467
468 Amg::Vector3D entryPoint(localStartPosition.x(),localStartPosition.y(),localStartPosition.z());
469 Amg::Vector3D exitPoint(localEndPosition.x(),localEndPosition.y(),localEndPosition.z());
470
471 const InDetDD::PixelModuleDesign* design(
dynamic_cast<const InDetDD::PixelModuleDesign*
>(&hitSiDetElement->
design()));
472 if (not design){
473 ATH_MSG_FATAL(
"Failed to cast the hitSiDetElement::design pointer to a PixelModuleDesign*; aborting due to null pointer");
474 return StatusCode::FAILURE;
475 }
476 const Amg::Vector2D localEntry(localStartPosition.x(),localStartPosition.y());
477 const Amg::Vector2D localExit(localEndPosition.x(),localEndPosition.y());
478
481
484
485 double halfthickness = hitSiDetElement->
thickness()*0.5;
486
487 bool EntryValid(entryCellId.
isValid());
488 bool ExitValid(exitCellId.
isValid());
489
491
492 Identifier diodeID = hitId;
495
497
498
499
500
501
502
503 pixMinimalPathCut=th0;
504
505 if (!EntryValid || !ExitValid)
506 {
507
508
509 if ( !EntryValid && !ExitValid) continue;
510
512
514
516 {
517 ATH_MSG_WARNING(
"ATTENTION THE INTERSECTION COORDINATE IS OUT OF THE MODULE");
518 continue;
519 }
520
521 const Amg::Vector2D Intersection_2d(Intersection.x(),Intersection.y());
523
524 InDetDD::SiCellId Intersection_2dCellId = hitSiDetElement->
cellIdFromIdentifier(Intersection_2dId);
525
526
527 if(!Intersection_2dCellId.
isValid())
continue;
528
529 if(EntryValid)
530 exitPoint = Intersection;
531 else
532 entryPoint = Intersection;
533
534
535 }
536
537
538
540 if(!digitizationModule){
542 return StatusCode::FAILURE;
543 }
544
545
546 std::vector<Trk::DigitizationStep> digitizationSteps =
m_digitizationStepper->cellSteps(*digitizationModule,entryPoint,exitPoint);
547
548
549
550
551 std::unique_ptr<InDet::PixelCluster>
pixelCluster =
nullptr;
553
554 std::vector<Identifier> rdoList;
555 std::vector<int> totList;
556
557 const bool isGanged = false;
558 int lvl1a = 0;
559
560 double accumulatedPathLength=0.;
561
562
563 int phiIndexMax = -999999;
564 int phiIndexMin = 1000000;
565 int etaIndexMax = -999999;
566 int etaIndexMin = 1000000;
567
568
569 for (auto& dStep : digitizationSteps){
570
572
574
579 }
580
581
582 if (pathlength < pixMinimalPathCut) continue;
583
584
587 InDetDD::SiCellId diode = hitSiDetElement->
cellIdOfPosition(PositionOnModule);
588
589
591 continue;
592
594
596 clusterPosition +=
pathlength * chargeCenterPosition;
597
598 int currentEtaIndex = diode.
etaIndex();
599 int currentPhiIndex = diode.
phiIndex();
600 if(currentEtaIndex > etaIndexMax) etaIndexMax = currentEtaIndex;
601 if(currentEtaIndex < etaIndexMin) etaIndexMin = currentEtaIndex;
602 if(currentPhiIndex > phiIndexMax) phiIndexMax = currentPhiIndex;
603 if(currentPhiIndex < phiIndexMin) phiIndexMin = currentPhiIndex;
604
605
606
608
609 rdoList.push_back(rdoId);
611
612 }
613
614 delete digitizationModule;
615
616
617 int siDeltaPhiCut = phiIndexMax-phiIndexMin+1;
618 int siDeltaEtaCut = etaIndexMax-etaIndexMin+1;
619
620 int totalToT=std::accumulate(totList.begin(), totList.end(), 0);;
621
622
623 if (rdoList.empty() || accumulatedPathLength < pixMinimalPathCut || totalToT == 0) {
624 if (totalToT == 0 && !rdoList.empty() )
ATH_MSG_WARNING(
"The total ToT of the cluster is 0, this should never happen");
625 continue;
626 }
627
628
629
630
631 clusterPosition *= 1./accumulatedPathLength;
633
634
635
636
637 bool merged = false;
639
640 for(Pixel_detElement_RIO_map::iterator currentClusIter = PixelDetElClusterMap.begin(); currentClusIter != PixelDetElClusterMap.end();) {
641
642 Pixel_detElement_RIO_map::iterator clusIter = currentClusIter++;
643 InDet::PixelCluster* currentCluster = clusIter->second;
644 const std::vector<Identifier> ¤tRdoList = currentCluster->
rdoList();
645 bool areNb = false;
646 for (auto rdoIter : rdoList) {
648 if (areNb) { break; }
649 }
650 if (areNb) {
651 const std::vector<int> ¤tTotList = currentCluster->
totList();
652 rdoList.insert(rdoList.end(), currentRdoList.begin(), currentRdoList.end() );
653 totList.insert(totList.end(), currentTotList.begin(), currentTotList.end() );
655 float c1 = (
float)currentRdoList.size();
656 float c2 = (
float)rdoList.size();
657 clusterPosition = (clusterPosition*
c2 + currentClusterPosition*
c1)/((c1 + c2));
659 merged = true;
660 PixelDetElClusterMap.erase(clusIter);
661
662
663 std::pair<PRD_MultiTruthCollection::iterator,PRD_MultiTruthCollection::iterator> saved_hit =
m_pixPrdTruth->equal_range(currentCluster->
identify());
664 for (PRD_MultiTruthCollection::iterator this_hit = saved_hit.first; this_hit != saved_hit.second; ++this_hit)
665 {
666 hit_vector.push_back(this_hit->second);
667 }
668
670 delete currentCluster;
671
672 }
673 }
674 }
675
676 bool not_valid = false;
677 for (auto & entry : rdoList) {
678 if (!(
entry.is_valid())) { not_valid =
true;
break;}
679 }
680
681 if (not_valid) continue;
682
683 if(merged) {
684
685 for (auto rdoIter : rdoList) {
687
688 int chargePhiIndex = chargeCellId.
phiIndex();
689 int chargeEtaIndex = chargeCellId.
etaIndex();
690
691 phiIndexMin = chargePhiIndex < phiIndexMin ? chargePhiIndex : phiIndexMin;
692 phiIndexMax = chargePhiIndex > phiIndexMax ? chargePhiIndex : phiIndexMax;
693 etaIndexMin = chargeEtaIndex < etaIndexMin ? chargeEtaIndex : etaIndexMin;
694 etaIndexMax = chargeEtaIndex > etaIndexMax ? chargeEtaIndex : etaIndexMax;
695 }
696 siDeltaPhiCut = (phiIndexMax-phiIndexMin)+1;
697 siDeltaEtaCut = (etaIndexMax-etaIndexMin)+1;
698 }
699
700
701
702
703
704
705
706
707
708
710
711
712 const auto* pixModDesign =
713 dynamic_cast<const InDetDD::PixelModuleDesign*>(
714 &hitSiDetElement->
design());
715 if (!pixModDesign) {
716 return StatusCode::FAILURE;
717 }
718 double etaWidth =
719 pixModDesign->widthFromColumnRange(etaIndexMin, etaIndexMax);
720
722 pixModDesign->widthFromRowRange(phiIndexMin, phiIndexMax);
723
724 InDet::SiWidth siWidth(
Amg::Vector2D(siDeltaPhiCut, siDeltaEtaCut),
726
727
730 clusterId, clusterPosition, std::move(rdoList),
731 lvl1a, std::move(totList), siWidth,
733 false, 0.0, 0.0, calibData, *offlineCalibData, ctx));
734 if (isGanged) {
736 }
737 } else {
739 "[ cluster - pix ] No pixels errors provided, but configured to "
740 "use them.");
742 " -> No pixels cluster will be created.");
743 continue;
744 }
745
748 continue;
749 }
750
753 continue;
754 }
755
757 PixelDetElClusterMap.insert(Pixel_detElement_RIO_map::value_type(
760
764 std::make_pair(insertedCluster->
identify(), currentLink));
766 << insertedCluster
767 << " and link = " << currentLink);
768 }
769 } else {
771 "Particle link NOT valid!! Truth map NOT filled with cluster"
772 << insertedCluster << " and link = " << currentLink);
773 }
774
775
776 for (const HepMcParticleLink& p : hit_vector) {
777
779 }
780
781 hit_vector.clear();
782 }
783
784
785 (void)
m_pixelClusterMap->insert(PixelDetElClusterMap.begin(), PixelDetElClusterMap.end());
786
787
788 }
789
790
791 return StatusCode::SUCCESS;
792}
#define ATH_MSG_ERROR(x,...)
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
int id() const
Return the id of the target particle.
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.
int barcode() const
Return the barcode of the target particle.
int phiIndex() const
Get phi index. Equivalent to strip().
bool isValid() const
Test if its in a valid state.
int etaIndex() const
Get eta index.
virtual SiCellId cellIdFromIdentifier(const Identifier &identifier) const override final
SiCellId from Identifier.
HepGeom::Point3D< double > hitLocalToLocal3D(const HepGeom::Point3D< double > &hitPosition) const
Same as previuos method but 3D.
SiCellId cellIdOfPosition(const Amg::Vector2D &localPos) const
As in previous method but returns SiCellId.
Identifier identifierOfPosition(const Amg::Vector2D &localPos) const
Full identifier of the cell for a given position: assumes a raw local position (no Lorentz shift).
Amg::Vector2D rawLocalPositionOfCell(const SiCellId &cellId) const
Returns position (center) of cell.
const std::vector< int > & totList() const
PixelChargeCalib::Thresholds getThresholds(InDetDD::PixelDiodeType type, unsigned int moduleHash, unsigned int FE) const
bool nextDetectorElement(const_iterator &b, const_iterator &e)
sets an iterator range with the hits of current detector element returns a bool when done
TimedVector::const_iterator const_iterator
const Segmentation & segmentation() const
return the segmenation
const Amg::Vector2D & localPosition() const
return the local position reference
Identifier identify() const
return the identifier
const std::vector< Identifier > & rdoList() const
return the List of rdo identifiers (pointers)
virtual const Amg::Vector2D cellPosition(const DigitizationCell &dCell) const =0
calculate the cell Position from the Id
Eigen::Matrix< double, 2, 1 > Vector2D
bool ignoreTruthLink(const T &p, bool vetoPileUp)
Helper function for SDO creation in PileUpTools.
std::pair< size_t, size_t > DigitizationCell
double pathlength() const
Method to retrieve the pathlength propagated till the Intersection.
PixelCluster_v1 PixelCluster
Define the version of the pixel cluster class.