Digitization functionality shared with RPC_PileUpTool.
Use special jitter consant for BIS & BIL chambers.
483 {
484 ATHRNG::RNGWrapper* rngWrapper =
m_rndmSvc->getEngine(
this);
486 CLHEP::HepRandomEngine* rndmEngine = rngWrapper->
getEngine(ctx);
487
488 const MuonGM::MuonDetectorManager* detMgr{nullptr};
490
491
492 std::unique_ptr<RPCSimHitCollection> inputSimHitColl{std::make_unique<RPCSimHitCollection>("RPC_Hits")};
493
494
495
496
498
499
502 return StatusCode::FAILURE;
503 }
504
505 struct SimDataContent {
507 std::vector<MuonSimData::Deposit> deposits;
509 double simTime{0.};
510 };
511
512 while (
m_thpcRPC->nextDetectorElement(i, e)) {
513
514
515 std::map<Identifier, SimDataContent> channelSimDataMap;
516
517
518 while (i != e) {
520
521 TimedHitPtr<RPCSimHit> phit(*i++);
522
523
524 const RPCSimHit& hit(*phit);
525
526 const int idHit = hit.RPCid();
527
528 const double globalHitTime{
hitTime(phit)};
529
530 const double G4Time{hit.globalTime()};
531
532 const double bunchTime{globalHitTime - hit.globalTime()};
533
534 ATH_MSG_DEBUG(
"Global time " << globalHitTime <<
" G4 time " << G4Time <<
" Bunch time " << bunchTime);
535
537 ATH_MSG_VERBOSE(
"Validation: globalHitTime, G4Time, BCtime = " << globalHitTime <<
" " << G4Time <<
" " << bunchTime);
538 inputSimHitColl->Emplace(idHit, globalHitTime, hit.localPosition(),
540 hit.postLocalPosition(),
541 hit.energyDeposit(), hit.stepLength(), hit.particleEncoding(), hit.kineticEnergy());
542 }
543
544
552
554
555
557 const Identifier elementID =
m_idHelper->elementID(stationName,stationEta,stationPhi,doubletR,
isValid);
559 ATH_MSG_WARNING(
"Failed to construct the element ID from "<<stationName
560 <<", stationEta: "<<stationEta<<", stationPhi: "<<stationPhi<<", doubletR: "<<doubletR);
561 continue;
562 }
563
565 << " stationName " << stationName << " stationEta " << stationEta << " stationPhi " << stationPhi << " doubletR "
566 << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi << " gasGap " << gasGap);
567 const Identifier detElId{
m_idHelper->channelID(elementID, doubletZ, doubletPhi, 1,0, 1,
isValid)};
569 continue;
570 }
575 }
576
577
578 bool isValidEta{false}, isValidPhi{false};
579 const Identifier idpaneleta =
m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 0, 1, isValidEta);
580 const Identifier idpanelphi =
m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 1, 1, isValidPhi);
581 if (!isValidEta || !isValidPhi) {
583 << " stationName " << stationName << " stationEta " << stationEta << " stationPhi " << stationPhi
584 << " doubletR " << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi << " gasGap "
585 << gasGap);
586 continue;
587 }
588
589
593 const double corrtimejitter = tmp_CorrJitter > 0.01 ?
594 CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0., tmp_CorrJitter) : 0.;
595
596
597
598
599
600 const Amg::Vector3D hitDir{(hit.postLocalPosition() - hit.localPosition()).unit()};
602 Amg::intersect<3>(hit.localPosition(), hitDir, Amg::Vector3D::UnitX(), 0).value_or(0) * hitDir;
603
604 std::array<int, 3> pcseta =
physicalClusterSize(ctx, reEle, idpaneleta, gapCentre, rndmEngine);
605 ATH_MSG_VERBOSE(
"Simulated cluster on eta panel: size/first/last= " << pcseta[0] <<
"/" << pcseta[1] <<
"/" << pcseta[2]);
606 std::array<int, 3> pcsphi =
physicalClusterSize(ctx, reEle, idpanelphi, gapCentre, rndmEngine);
607 ATH_MSG_VERBOSE(
"Simulated cluster on phi panel: size/first/last= " << pcsphi[0] <<
"/" << pcsphi[1] <<
"/" << pcsphi[2]);
608
609
610
611
612 const Identifier atlasRpcIdeta =
m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 0, pcseta[1], isValidEta);
613 const Identifier atlasRpcIdphi =
m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 1, pcsphi[1], isValidPhi);
614
616 const auto [etaStripOn, phiStripOn] =
detectionEfficiency(ctx, idpaneleta, idpanelphi, rndmEngine, particleLink);
617 ATH_MSG_DEBUG(
"SetPhiOn " << phiStripOn <<
" SetEtaOn " << etaStripOn);
618
619 for (bool imeasphi : {false, true}) {
620 if (!imeasphi && (!etaStripOn || !isValidEta)) continue;
621 if (imeasphi && (!phiStripOn || !isValidPhi)) continue;
622
623
624
625 const Identifier& atlasId = !imeasphi ? atlasRpcIdeta : atlasRpcIdphi;
626 std::array<int, 3>
pcs{!imeasphi ? pcseta : pcsphi};
627
628 ATH_MSG_DEBUG(
"SetOn: stationName " << stationName <<
" stationEta " << stationEta <<
" stationPhi " << stationPhi
629 << " doubletR " << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi
630 << " gasGap " << gasGap << " measphi " << imeasphi);
631
632
634 if (pcs[2] < 0){
635 continue;
636 }
637
638 ATH_MSG_DEBUG(
"Simulated cluster1: size/first/last= " << pcs[0] <<
"/" << pcs[1] <<
"/" << pcs[2]);
639
640
643
645 <<
" hit "<<
m_idHelper->print_to_string(atlasId)
653
654
656
657 double tns = G4Time + proptime + corrtimejitter;
658 ATH_MSG_VERBOSE(
"TOF+propagation time " << tns <<
" /s where proptime " << proptime <<
"/s");
659
660 double time = tns + bunchTime;
662
663
665
666 double*
b =
reinterpret_cast<double*
>(&packedMCword);
667
669
670
671
673
674
675
677 if (std::abs(hit.particleEncoding()) == 13 || hit.particleEncoding() == 0) {
678 if (channelSimDataMap.find(atlasId) == channelSimDataMap.end()) {
679 SimDataContent&
content = channelSimDataMap[atlasId];
681 content.deposits.push_back(deposit);
686 }
687 }
688 }
689
690
691
692
693
694
695
697 for (
int clus = pcs[1]; clus <=
pcs[2]; ++clus) {
698 Identifier newId =
m_idHelper->channelID(stationName, stationEta, stationPhi, doubletR, doubletZ,
699 doubletPhi, gasGap, imeasphi, clus,
isValid);
701 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
"Channel "<< stationName<<
" "<<stationEta<<
" "<<stationPhi<<
" "<< doubletR<<
" "<<doubletZ
702 <<" "<< doubletPhi<<" "<< gasGap <<" "<< imeasphi<<" "<< clus<<" is invalid");
703 continue;
704 }
705
708 ATH_MSG_WARNING(
"Temporary skipping creation of RPC digit for stationName="
709 << stationName << ", eta=" << stationEta << ", phi=" << stationPhi << ", doubletR=" << doubletR
710 << ", doubletZ=" << doubletZ << ", doubletPhi=" << doubletPhi << ", gasGap=" << gasGap
711 << ", measuresPhi=" << imeasphi << ", strip=" << clus << ", cf. ATLASRECTS-6124");
712 return StatusCode::SUCCESS;
713 } else {
716 return StatusCode::FAILURE;
717 }
718 }
719
723
725 std::vector<MuonSimData::Deposit> newdeps;
726 newdeps.push_back(deposit);
727 m_sdo_tmp_map.insert(std::map<Identifier, std::vector<MuonSimData::Deposit>>::value_type(newId, newdeps));
728 } else {
730 }
731 }
732 }
733 }
734
736 for (auto it = channelSimDataMap.begin(); it != channelSimDataMap.end(); ++it) {
737 MuonSimData
simData(
it->second.deposits, 0);
740 auto insertResult = sdoContainer->insert(std::make_pair(
it->first,
simData));
741 if (!insertResult.second)
742 ATH_MSG_WARNING(
"Attention: this sdo is not recorded, since the identifier already exists in the sdoContainer map");
743 }
744 }
745
746 }
747
749
750 std::map<Identifier, std::vector<MuonSimData::Deposit>>::iterator map_iter =
m_sdo_tmp_map.begin();
752
754
755 const Identifier theId = (*map_iter).first;
757
758 const std::vector<MuonSimData::Deposit> theDeps = (*map_iter).second;
759
760
762 std::multimap<double, MuonSimData::Deposit>
times;
763
764
765 for (
unsigned int k = 0;
k < theDeps.size();
k++) {
766 double time = theDeps[
k].second.secondEntry();
767 times.insert(std::multimap<double, MuonSimData::Deposit>::value_type(time, theDeps[k]));
768 }
769
770
771
772 IdContext rpcContext =
m_idHelper->module_context();
773
774 std::multimap<double, MuonSimData::Deposit>::iterator map_dep_iter =
times.begin();
775
776
777 double last_time = -10000;
778 for (; map_dep_iter !=
times.end(); ++map_dep_iter) {
779 double currTime = (*map_dep_iter).first;
781
783
784
785 if (sdoContainer->find(theId) != sdoContainer->end())
786 {
787 std::map<Identifier, MuonSimData>::const_iterator
it = sdoContainer->find(theId);
788 std::vector<MuonSimData::Deposit> deps = ((*it).second).getdeposits();
789 deps.push_back((*map_dep_iter).second);
790 } else
791 {
792 std::vector<MuonSimData::Deposit> deposits;
793 deposits.push_back((*map_dep_iter).second);
794 std::pair<std::map<Identifier, MuonSimData>::iterator, bool> insertResult =
795 sdoContainer->insert(std::make_pair(theId, MuonSimData(deposits, 0)));
796 if (!insertResult.second)
798 "Attention TEMP: this sdo is not recorded, since the identifier already exists in the sdoContainer map");
799 }
800 }
801
802 if (std::abs(currTime - last_time) > (
m_deadTime)) {
803 ATH_MSG_DEBUG(
"deposit with time " << currTime <<
" is distant enough from previous (if any) hit on teh same strip");
804 last_time = (*map_dep_iter).first;
805
806
807 double uncorrjitter = 0;
810 if (tmp_UncorrJitter > 0.01) uncorrjitter = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0., tmp_UncorrJitter);
811
812
813
814
818
820 double newDigit_time = currTime + uncorrjitter +
m_rpc_time_shift -
tp - propTimeFromStripCenter;
821
822 double digi_ToT = -1.;
824
825 ATH_MSG_VERBOSE(
"last_time=currTime " << last_time <<
" jitter " << uncorrjitter <<
" TOFcorrection " << tp <<
" shift "
827
828
829 bool outsideDigitizationWindow =
outsideWindow(newDigit_time);
830 if (outsideDigitizationWindow) {
831 ATH_MSG_VERBOSE(
"hit outside digitization window - do not produce digits");
835
836 continue;
837 }
838
839
840 last_time = (*map_dep_iter).first;
841
842 std::unique_ptr<RpcDigit> newDigit = std::make_unique<RpcDigit>(theId, newDigit_time, digi_ToT, false);
843
844 Identifier elemId =
m_idHelper->elementID(theId);
845 RpcDigitCollection* digitCollection = nullptr;
846
847 IdentifierHash coll_hash;
848 if (
m_idHelper->get_hash(elemId, coll_hash, &rpcContext)) {
849 ATH_MSG_ERROR(
"Unable to get RPC hash id from RPC Digit collection "
850 <<
"context begin_index = " << rpcContext.
begin_index()
851 <<
" context end_index = " << rpcContext.
end_index() <<
" the identifier is ");
853 }
854
855
857
858
859 if (coll_hash >= collections.size()) {
860 collections.resize (coll_hash+1);
861 }
862 digitCollection = collections[coll_hash].
get();
863 if (!digitCollection) {
864 collections[coll_hash] = std::make_unique<RpcDigitCollection>(elemId, coll_hash);
865 digitCollection = collections[coll_hash].
get();
866 }
867 digitCollection->
push_back(std::move(newDigit));
868
870
871 if (sdoContainer->find(theId) != sdoContainer->end()) {
872 std::map<Identifier, MuonSimData>::const_iterator
it = sdoContainer->find(theId);
873 std::vector<MuonSimData::Deposit> deps = ((*it).second).getdeposits();
874 deps.push_back((*map_dep_iter).second);
875 } else {
876 std::vector<MuonSimData::Deposit> deposits;
877 deposits.push_back((*map_dep_iter).second);
878 std::pair<std::map<Identifier, MuonSimData>::iterator, bool> insertResult =
879 sdoContainer->insert(std::make_pair(theId, MuonSimData(deposits, 0)));
880 if (!insertResult.second)
882 "Attention: this sdo is not recorded, since teh identifier already exists in the sdoContainer map");
883 }
884 }
885
886 } else
887 ATH_MSG_DEBUG(
"discarding digit due to dead time: " << (*map_dep_iter).first <<
" " << last_time);
888 }
889
890 }
891
892
897 }
898
899 return StatusCode::SUCCESS;
900}
float hitTime(const AFP_SIDSimHit &hit)
#define ATH_CHECK
Evaluate an expression and check for errors.
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
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.
const T * get(size_type n) const
Access an element, as an rvalue.
value_type push_back(value_type pElem)
Add an element to the end of the collection.
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.
size_type begin_index() const
size_type end_index() const
void show() const
Print out in hex form.
virtual const Amg::Transform3D & transform() const override
Return local to global transform.
const RpcReadoutElement * getRpcReadoutElement(const Identifier &id) const
access via extended identifier (requires unpacking)
bool rotatedRpcModule() const
Amg::Vector3D stripPos(const Identifier &id) const
Amg::Transform3D localToGlobalTransf(const Identifier &id) const
std::pair< HepMcParticleLink, MuonMCData > Deposit
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
TimedVector::const_iterator const_iterator
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Eigen::Matrix< double, 3, 1 > Vector3D
time(flags, cells_name, *args, **kw)
bool ignoreTruthLink(const T &p, bool vetoPileUp)
Helper function for SDO creation in PileUpTools.
constexpr uint8_t stationPhi
station Phi 1 to 8