405{
406
409 rngWrapper->
setSeed( rngName, ctx );
410 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
411
412 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
m_pixelDetEleCollKey, ctx);
413 const InDetDD::SiDetectorElementCollection* elements(*pixelDetEleHandle);
414 if (not pixelDetEleHandle.isValid() or elements==nullptr) {
416 return StatusCode::FAILURE;
417 }
418
420
423
424 SG::ReadCondHandle<PixelChargeCalibCondData> calibDataHandle(
m_chargeDataKey, ctx);
425 const PixelChargeCalibCondData *calibData = *calibDataHandle;
427 std::vector<int> truthIdList;
428 std::vector<Identifier> detEl;
429
431
433
434 truthIdList.clear();
435 detEl.clear();
436
437 while (i != e) {
438
439 const TimedHitPtr<SiHit>& hit(*i++);
440
441
443 const int barrelEC = hit->getBarrelEndcap();
444 const int layerDisk = hit->getLayerDisk();
445 const int phiModule = hit->getPhiModule();
446 const int etaModule = hit->getEtaModule();
447
448 const Identifier moduleID =
m_pixel_ID->wafer_id(barrelEC, layerDisk, phiModule, etaModule);
449 const IdentifierHash waferHash =
m_pixel_ID->wafer_hash(moduleID);
450 const InDetDD::SiDetectorElement* hitSiDetElement = elements->getDetectorElement(waferHash);
451 if (!hitSiDetElement) {
ATH_MSG_ERROR(
" could not get detector element ");
continue;}
452
453 if (!(hitSiDetElement->
isPixel())) {
continue;}
454
455
456
457 std::vector<HepMcParticleLink> hit_vector;
458
459 const int truthID = (currentLink.
barcode() !=0 && currentLink.
id() == 0) ? 3 : currentLink.
id();
460
461 const Identifier hitId = hitSiDetElement->
identify();
462
463
464
465 bool isRep = false;
466
467 for (int j : truthIdList) {
468 for (auto & k : detEl) {
469 if ((truthID > 0) && (truthID == j) && (hitId == k)) {isRep = true; break;}
470 }
471 if (isRep) break;
472 }
473
474 if (isRep) continue;
475
476 truthIdList.push_back(truthID);
477 detEl.push_back(hitId);
478
479 HepGeom::Point3D<double> localStartPosition = hit->localStartPosition();
480 HepGeom::Point3D<double> localEndPosition = hit->localEndPosition();
481
484
485 int isEndcap = (barrelEC==0) ? 0:1;
486
489
490
491
492
493 Diffuse(localStartPosition, localEndPosition, shiftX, shiftY);
494
495 const Amg::Vector3D localDirection(localEndPosition.x()-localStartPosition.x(), localEndPosition.y()-localStartPosition.y(), localEndPosition.z()-localStartPosition.z());
496
497 Amg::Vector3D entryPoint(localStartPosition.x(),localStartPosition.y(),localStartPosition.z());
498 Amg::Vector3D exitPoint(localEndPosition.x(),localEndPosition.y(),localEndPosition.z());
499
500 const InDetDD::PixelModuleDesign* design(
dynamic_cast<const InDetDD::PixelModuleDesign*
>(&hitSiDetElement->
design()));
501 if (not design){
502 ATH_MSG_FATAL(
"Failed to cast the hitSiDetElement::design pointer to a PixelModuleDesign*; aborting due to null pointer");
503 return StatusCode::FAILURE;
504 }
505 const Amg::Vector2D localEntry(localStartPosition.x(),localStartPosition.y());
506 const Amg::Vector2D localExit(localEndPosition.x(),localEndPosition.y());
507
510
513
514 double halfthickness = hitSiDetElement->
thickness()*0.5;
515
516 bool EntryValid(entryCellId.
isValid());
517 bool ExitValid(exitCellId.
isValid());
518
520
521 Identifier diodeID = hitId;
524
526
527
528
529
530
531
532 pixMinimalPathCut=th0;
533
534 if (!EntryValid || !ExitValid)
535 {
536
537
538 if ( !EntryValid && !ExitValid) continue;
539
541
543
545 {
546 ATH_MSG_WARNING(
"ATTENTION THE INTERSECTION COORDINATE IS OUT OF THE MODULE");
547 continue;
548 }
549
550 const Amg::Vector2D Intersection_2d(Intersection.x(),Intersection.y());
552
553 InDetDD::SiCellId Intersection_2dCellId = hitSiDetElement->
cellIdFromIdentifier(Intersection_2dId);
554
555
556 if(!Intersection_2dCellId.
isValid())
continue;
557
558 if(EntryValid)
559 exitPoint = Intersection;
560 else
561 entryPoint = Intersection;
562
563
564 }
565
566
567
569 if(!digitizationModule){
571 return StatusCode::FAILURE;
572 }
573
574
575 std::vector<Trk::DigitizationStep> digitizationSteps =
m_digitizationStepper->cellSteps(*digitizationModule,entryPoint,exitPoint);
576
577
578
579
580 std::unique_ptr<InDet::PixelCluster>
pixelCluster =
nullptr;
582
583 std::vector<Identifier> rdoList;
584 std::vector<int> totList;
585
586 const bool isGanged = false;
587 int lvl1a = 0;
588
589 double accumulatedPathLength=0.;
590
591
592 int phiIndexMax = -999999;
593 int phiIndexMin = 1000000;
594 int etaIndexMax = -999999;
595 int etaIndexMin = 1000000;
596
597
598 for (auto& dStep : digitizationSteps){
599
600 double pathlength = dStep.stepLength;
601
603
607 pathlength *= (1.+sPar);
608 }
609
610
611 if (pathlength < pixMinimalPathCut) continue;
612
613
616 InDetDD::SiCellId diode = hitSiDetElement->
cellIdOfPosition(PositionOnModule);
617
618
620 continue;
621
623
625 clusterPosition += pathlength * chargeCenterPosition;
626
627 int currentEtaIndex = diode.
etaIndex();
628 int currentPhiIndex = diode.
phiIndex();
629 if(currentEtaIndex > etaIndexMax) etaIndexMax = currentEtaIndex;
630 if(currentEtaIndex < etaIndexMin) etaIndexMin = currentEtaIndex;
631 if(currentPhiIndex > phiIndexMax) phiIndexMax = currentPhiIndex;
632 if(currentPhiIndex < phiIndexMin) phiIndexMin = currentPhiIndex;
633
634
635
636 accumulatedPathLength += pathlength;
637
638 rdoList.push_back(rdoId);
640
641 }
642
643 delete digitizationModule;
644
645
646 int siDeltaPhiCut = phiIndexMax-phiIndexMin+1;
647 int siDeltaEtaCut = etaIndexMax-etaIndexMin+1;
648
649 int totalToT=std::accumulate(totList.begin(), totList.end(), 0);;
650
651
652 if (rdoList.empty() || accumulatedPathLength < pixMinimalPathCut || totalToT == 0) {
653 if (totalToT == 0 && !rdoList.empty() )
ATH_MSG_WARNING(
"The total ToT of the cluster is 0, this should never happen");
654 continue;
655 }
656
657
658 clusterPosition *= 1./accumulatedPathLength;
660
661
662
663
664 bool merged = false;
666
667 for(Pixel_detElement_RIO_map::iterator currentClusIter = PixelDetElClusterMap.begin(); currentClusIter != PixelDetElClusterMap.end();) {
668
669 Pixel_detElement_RIO_map::iterator clusIter = currentClusIter++;
670 InDet::PixelCluster* currentCluster = clusIter->second;
671 const std::vector<Identifier> ¤tRdoList = currentCluster->
rdoList();
672 bool areNb = false;
673 for (auto rdoIter : rdoList) {
675 if (areNb) { break; }
676 }
677 if (areNb) {
678 const std::vector<int> ¤tTotList = currentCluster->
totList();
679 rdoList.insert(rdoList.end(), currentRdoList.begin(), currentRdoList.end() );
680 totList.insert(totList.end(), currentTotList.begin(), currentTotList.end() );
682 float c1 = (
float)currentRdoList.size();
683 float c2 = (
float)rdoList.size();
684 clusterPosition = (clusterPosition*
c2 + currentClusterPosition*
c1)/((c1 + c2));
686 merged = true;
687 PixelDetElClusterMap.erase(clusIter);
688
689
690 std::pair<PRD_MultiTruthCollection::iterator,PRD_MultiTruthCollection::iterator> saved_hit =
m_pixPrdTruth->equal_range(currentCluster->
identify());
691 for (PRD_MultiTruthCollection::iterator this_hit = saved_hit.first; this_hit != saved_hit.second; ++this_hit)
692 {
693 hit_vector.push_back(this_hit->second);
694 }
695
697 delete currentCluster;
698
699 }
700 }
701 }
702
703 bool not_valid = false;
704 for (auto & entry : rdoList) {
705 if (!(
entry.is_valid())) { not_valid =
true;
break;}
706 }
707
708 if (not_valid) continue;
709
710 if(merged) {
711
712 for (auto rdoIter : rdoList) {
714
715 int chargePhiIndex = chargeCellId.
phiIndex();
716 int chargeEtaIndex = chargeCellId.
etaIndex();
717
718 phiIndexMin = chargePhiIndex < phiIndexMin ? chargePhiIndex : phiIndexMin;
719 phiIndexMax = chargePhiIndex > phiIndexMax ? chargePhiIndex : phiIndexMax;
720 etaIndexMin = chargeEtaIndex < etaIndexMin ? chargeEtaIndex : etaIndexMin;
721 etaIndexMax = chargeEtaIndex > etaIndexMax ? chargeEtaIndex : etaIndexMax;
722 }
723 siDeltaPhiCut = (phiIndexMax-phiIndexMin)+1;
724 siDeltaEtaCut = (etaIndexMax-etaIndexMin)+1;
725 }
726
727
728
729
730
731
732
733
734
735
737
738
739 const auto* pixModDesign =
740 dynamic_cast<const InDetDD::PixelModuleDesign*>(
741 &hitSiDetElement->
design());
742 if (!pixModDesign) {
743 return StatusCode::FAILURE;
744 }
745 double etaWidth =
746 pixModDesign->widthFromColumnRange(etaIndexMin, etaIndexMax);
747
749 pixModDesign->widthFromRowRange(phiIndexMin, phiIndexMax);
750
751 InDet::SiWidth siWidth(
Amg::Vector2D(siDeltaPhiCut, siDeltaEtaCut),
753
754
757 clusterId, clusterPosition, std::vector<Identifier>(rdoList),
758 lvl1a, std::vector<int>(totList), siWidth,
760 false, 0.0, 0.0, calibData, *offlineCalibData, ctx));
761 if (isGanged) {
763 }
764 } else {
766 "[ cluster - pix ] No pixels errors provided, but configured to "
767 "use them.");
769 " -> No pixels cluster will be created.");
770 continue;
771 }
772
775 continue;
776 }
777
780 continue;
781 }
782
784 PixelDetElClusterMap.insert(Pixel_detElement_RIO_map::value_type(
787
791 std::make_pair(insertedCluster->
identify(), currentLink));
793 << insertedCluster
794 << " and link = " << currentLink);
795 }
796 } else {
798 "Particle link NOT valid!! Truth map NOT filled with cluster"
799 << insertedCluster << " and link = " << currentLink);
800 }
801
802
803 for (const HepMcParticleLink& p : hit_vector) {
804
806 }
807
808 hit_vector.clear();
809 }
810
811
812 (void)
m_pixelClusterMap->insert(PixelDetElClusterMap.begin(), PixelDetElClusterMap.end());
813
814
815 }
816
817
818 return StatusCode::SUCCESS;
819}
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
PixelCluster_v1 PixelCluster
Define the version of the pixel cluster class.