diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonTruthAlgsR4/python/MuonTruthAlgsConfig.py b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonTruthAlgsR4/python/MuonTruthAlgsConfig.py
index 94aa98d122f58952acf5a781f23987a889941ecb..2dfee864ff8874e42d47093303ce396ab566db65 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonTruthAlgsR4/python/MuonTruthAlgsConfig.py
+++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonTruthAlgsR4/python/MuonTruthAlgsConfig.py
@@ -17,6 +17,12 @@ def TruthSegmentMakerCfg(flags, name = "TruthSegmentMakerAlg", **kwargs):
if flags.Detector.EnableMM: containerNames+=["xMmSimHits"]
if flags.Detector.EnablesTGC: containerNames+=["xStgcSimHits"]
kwargs.setdefault("SimHitKeys", containerNames)
+ PrdLinkInputs = []
+ if flags.Detector.EnableMDT: PrdLinkInputs+=[( 'xAOD::UncalibratedMeasurementContainer' , 'StoreGateSvc+xAODMdtCircles.simHitLink' )]
+ if flags.Detector.EnableRPC: PrdLinkInputs+=[ ( 'xAOD::UncalibratedMeasurementContainer' , 'StoreGateSvc+xRpcMeasurements.simHitLink' )]
+ if flags.Detector.EnableTGC: PrdLinkInputs+=[('xAOD::UncalibratedMeasurementContainer' , 'StoreGateSvc+xTgcStrips.simHitLink' )]
+
+ kwargs.setdefault("ExtraInputs", PrdLinkInputs)
the_alg = CompFactory.MuonR4.TruthSegmentMaker(name, **kwargs)
result.addEventAlgo(the_alg, primary = True)
return result
diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.cxx b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.cxx
index df4ea07d2abcea06a8177216a2c7248030fc52dc..490bc7df5194dba0f13a624ab5a531f1ced92ecd 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.cxx
@@ -22,6 +22,8 @@ namespace MuonValR4 {
ATH_CHECK(m_inSimHitKeys.initialize());
ATH_CHECK(m_spacePointKey.initialize());
ATH_CHECK(m_inHoughSegmentSeedKey.initialize());
+ ATH_CHECK(m_truthSegmentKey.initialize());
+ ATH_CHECK(m_rh_truthSegmentSimHitLink.initialize());
ATH_CHECK(m_inSegmentKey.initialize(!m_inSegmentKey.empty()));
m_tree.addBranch(std::make_shared<MuonVal::EventInfoBranch>(m_tree,0));
m_out_SP = std::make_shared<MuonValR4::SpacePointTesterModule>(m_tree, m_spacePointKey.key());
@@ -64,30 +66,33 @@ namespace MuonValR4 {
}
- void MuonHoughTransformTester::matchSeedToTruth(const MuonR4::SegmentSeed* seed,
+ void MuonHoughTransformTester::matchSeedToTruth(const EventContext & ctx, const MuonR4::SegmentSeed* seed,
chamberLevelObjects & objs ) const{
double bestTruthFrac{0.};
- HepMC::ConstGenParticlePtr bestMatch = nullptr;
- for (auto & [ genParticle, truthQuantities] : objs.truthMatching) {
+ const xAOD::MuonSegment* bestMatch = nullptr;
+ SG::ReadDecorHandle<xAOD::MuonSegmentContainer,SimHitLinkVec> simHitsFromTruth(m_rh_truthSegmentSimHitLink, ctx);
+ SG::ConstAccessor<ElementLink<xAOD::MuonSimHitContainer>> simHitAcc("simHitLink");
+
+ for (auto & [ truthsegment, truthQuantities] : objs.truthMatching) {
unsigned int nRecFound{0};
+ const SimHitLinkVec & simHitsForSegment = simHitsFromTruth(*truthsegment);
+
for (const MuonR4::HoughHitType& spacePoint : seed->getHitsInMax()) {
- for (const xAOD::MuonSimHit* simHit : truthQuantities.detectorHits) {
- if(m_idHelperSvc->isMdt(simHit->identify()) &&
- spacePoint->identify() == simHit->identify()){
- ++nRecFound;
- break;
- }
- //// Sim hits are expressed w.r.t to the gas gap Id. Check whether
- /// the hit is in the same gas gap
- else if (m_idHelperSvc->gasGapId(spacePoint->identify()) == simHit->identify()) {
- ++nRecFound;
- // break; // should we... ?
- }
+ if (!simHitAcc.isAvailable(*spacePoint->primaryMeasurement())) continue;
+ auto simHitMatch = simHitAcc(*spacePoint->primaryMeasurement());
+ if (!simHitMatch.isValid() || *simHitMatch == nullptr){
+ continue;
+ }
+ auto found = std::ranges::find_if( simHitsForSegment, [&simHitMatch](const ElementLink<xAOD::MuonSimHitContainer> & el){
+ return (el == simHitMatch);
+ });
+ if (found != simHitsForSegment.end()){
+ ++nRecFound;
}
}
double truthFraction = (1.*nRecFound) / (1.*seed->getHitsInMax().size());
if (truthFraction > bestTruthFrac) {
- bestMatch = genParticle;
+ bestMatch = truthsegment;
bestTruthFrac = truthFraction;
}
}
@@ -95,14 +100,14 @@ namespace MuonValR4 {
/** Map the seed to the truth particle */
chamberLevelObjects::SeedMatchQuantites& seedMatch = objs.seedMatching[seed];
seedMatch.matchProb = bestTruthFrac;
- seedMatch.truthParticle = bestMatch;
+ seedMatch.truthsegment = bestMatch;
/** Back mapping of the best truth -> seed */
objs.truthMatching[bestMatch].assocSeeds.push_back(seed);
}
- void MuonHoughTransformTester::matchSeedsToTruth(chamberLevelObjects & objs) const {
+ void MuonHoughTransformTester::matchSeedsToTruth(const EventContext & ctx,chamberLevelObjects & objs) const {
for (auto & [ seed, matchObj] : objs.seedMatching) {
- matchSeedToTruth(seed, objs);
+ matchSeedToTruth(ctx, seed, objs);
ATH_MSG_VERBOSE("Truth matching probability "<<matchObj.matchProb);
}
}
@@ -112,37 +117,42 @@ namespace MuonValR4 {
m_out_stationEta = chamber->stationEta();
m_out_stationPhi = chamber->stationPhi();
}
- void MuonHoughTransformTester::fillTruthInfo(const HepMC::ConstGenParticlePtr genParticlePtr, const std::vector<const xAOD::MuonSimHit*> & hits,const ActsGeometryContext & gctx){
- if (!genParticlePtr) return;
+ void MuonHoughTransformTester::fillTruthInfo(const EventContext & ctx, const xAOD::MuonSegment* segment,const ActsGeometryContext & gctx){
+ if (!segment) return;
m_out_hasTruth = true;
- m_out_gen_Eta = genParticlePtr->momentum().eta();
- m_out_gen_Phi= genParticlePtr->momentum().phi();
- m_out_gen_Pt= genParticlePtr->momentum().perp();
-
- const xAOD::MuonSimHit* simHit = hits.front();
- const Identifier ID = simHit->identify();
-
- const Amg::Transform3D toChamber{toChamberTrf(gctx, ID)};
- const Amg::Vector3D localPos{toChamber * xAOD::toEigen(simHit->localPosition())};
- Amg::Vector3D chamberDir = toChamber.linear() * xAOD::toEigen(simHit->localDirection());
+ Amg::Vector3D segPos{
+ segment->x(),
+ segment->y(),
+ segment->z(),
+ };
+ Amg::Vector3D segDir{
+ segment->px(),
+ segment->py(),
+ segment->pz(),
+ };
+ // eta is interpreted as the eta-location
+ m_out_gen_Eta = segDir.eta();
+ m_out_gen_Phi= segDir.phi();
+ m_out_gen_Pt= segDir.perp();
+
+ SG::ReadDecorHandle<xAOD::MuonSegmentContainer,SimHitLinkVec> simHitsFromTruth(m_rh_truthSegmentSimHitLink, ctx);
+ const SimHitLinkVec & simHits = simHitsFromTruth(*segment);
+ const xAOD::MuonSimHit* firstSimHit = *(simHits.front());
+ const Identifier ID = firstSimHit->identify();
+ const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(ID);
+ //transform from local (w.r.t tube's frame) to global (ATLAS frame) and then to chamber's frame
+ const MuonGMR4::MuonChamber* muonChamber = reElement->getChamber();
+ auto toChamber = muonChamber->globalToLocalTrans(gctx);
+ const Amg::Vector3D chamberPos{toChamber * segPos};
+ Amg::Vector3D chamberDir = toChamber.linear() * segDir;
- /// Express the simulated hit in the center of the chamber
- const std::optional<double> lambda = Amg::intersect<3>(localPos, chamberDir, Amg::Vector3D::UnitZ(), 0.);
- Amg::Vector3D chamberPos = localPos + (*lambda)*chamberDir;
- m_out_gen_nHits = hits.size();
- unsigned int nMdt{0}, nRpc{0}, nTgc{0}, nMm{0}, nsTgc{0};
- for (const xAOD::MuonSimHit* hit : hits){
- nMdt += m_idHelperSvc->isMdt(hit->identify());
- nRpc += m_idHelperSvc->isRpc(hit->identify());
- nTgc += m_idHelperSvc->isTgc(hit->identify());
- nMm += m_idHelperSvc->isMM(hit->identify());
- nsTgc += m_idHelperSvc->issTgc(hit->identify());
- }
- m_out_gen_nRPCHits = nRpc;
- m_out_gen_nMDTHits = nMdt;
- m_out_gen_nTGCHits = nTgc;
- m_out_gen_nMMits = nMm;
- m_out_gen_nsTGCHits = nsTgc;
+ m_out_gen_nHits = segment->nPrecisionHits()+segment->nPhiLayers() + segment->nTrigEtaLayers();
+
+ m_out_gen_nMDTHits = (segment->technology() == Muon::MuonStationIndex::MDT ? segment->nPrecisionHits() : 0);
+ m_out_gen_nNswHits = (segment->technology() != Muon::MuonStationIndex::MDT ? segment->nPrecisionHits() : 0);
+ m_out_gen_nTGCHits = (segment->chamberIndex() > Muon::MuonStationIndex::ChIndex::BEE ? segment->nPhiLayers() + segment->nTrigEtaLayers() : 0);
+ m_out_gen_nRPCHits = (segment->chamberIndex() <= Muon::MuonStationIndex::ChIndex::BEE ? segment->nPhiLayers() + segment->nTrigEtaLayers() : 0);
+
m_out_gen_tantheta = (std::abs(chamberDir.z()) > 1.e-8 ? chamberDir.y()/chamberDir.z() : 1.e10);
m_out_gen_tanphi = (std::abs(chamberDir.z()) > 1.e-8 ? chamberDir.x()/chamberDir.z() : 1.e10);
@@ -239,6 +249,11 @@ namespace MuonValR4 {
const MuonR4::SegmentContainer* readMuonSegments{nullptr};
ATH_CHECK(retrieveContainer(ctx, m_inSegmentKey, readMuonSegments));
+
+ const xAOD::MuonSegmentContainer* readTruthSegments{nullptr};
+ ATH_CHECK(retrieveContainer(ctx, m_truthSegmentKey, readTruthSegments));
+
+ SG::ReadDecorHandle<xAOD::MuonSegmentContainer, SimHitLinkVec> readSimHits{m_rh_truthSegmentSimHitLink, ctx};
ATH_MSG_DEBUG("Succesfully retrieved input collections");
@@ -246,26 +261,18 @@ namespace MuonValR4 {
// The fast digi should only generate one circle per tube.
std::map<const MuonGMR4::MuonChamber*, chamberLevelObjects> allObjectsPerChamber;
- for (const SG::ReadHandleKey<xAOD::MuonSimHitContainer>& key : m_inSimHitKeys){
- const xAOD::MuonSimHitContainer* collection{nullptr};
- ATH_CHECK(retrieveContainer(ctx, key, collection));
- for (const xAOD::MuonSimHit* simHit : *collection) {
- const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(simHit->identify());
- const MuonGMR4::MuonChamber* id{reElement->getChamber()};
- chamberLevelObjects & theObjects = allObjectsPerChamber[id];
- auto genLink = simHit->genParticleLink();
- HepMC::ConstGenParticlePtr genParticle = nullptr;
- if (genLink.isValid()){
- genParticle = genLink.cptr();
- }
- /// skip empty truth matches for now
- if (!genParticle) continue;
- theObjects.truthMatching[genParticle].detectorHits.push_back(simHit);
- }
+ for (const xAOD::MuonSegment* truthSegment : *readTruthSegments){
+ const SimHitLinkVec & simHits = readSimHits(*truthSegment);
+ const Identifier simHitId = (*simHits.front())->identify();
+ const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(simHitId);
+ const MuonGMR4::MuonChamber* id{reElement->getChamber()};
+ chamberLevelObjects & theObjects = allObjectsPerChamber[id];
+ theObjects.truthMatching[truthSegment].truthsegment = truthSegment;
}
// Populate the seeds first
for (const MuonR4::SegmentSeed* max : *readSegmentSeeds) {
+ // use adventurous minion syntax to add a key with a default value
allObjectsPerChamber[max->chamber()].seedMatching[max];
}
if (readMuonSegments) {
@@ -277,18 +284,18 @@ namespace MuonValR4 {
}
for (auto & [chamber, chamberLevelObjects] : allObjectsPerChamber){
- matchSeedsToTruth(chamberLevelObjects);
+ matchSeedsToTruth(ctx, chamberLevelObjects);
/// Step 1: Fill the matched pairs
- for (auto & [genParticlePtr, assocInfo] : chamberLevelObjects.truthMatching) {
+ for (auto & [truthSegment, assocInfo] : chamberLevelObjects.truthMatching) {
if (assocInfo.assocSeeds.empty()) {
fillChamberInfo(chamber);
- fillTruthInfo(genParticlePtr, assocInfo.detectorHits, gctx);
+ fillTruthInfo(ctx, truthSegment, gctx);
if (!m_tree.fill(ctx)) return StatusCode::FAILURE;
continue;
}
for (const MuonR4::SegmentSeed* seed : assocInfo.assocSeeds) {
fillChamberInfo(chamber);
- fillTruthInfo(genParticlePtr, assocInfo.detectorHits, gctx);
+ fillTruthInfo(ctx, truthSegment, gctx);
auto& seedMatch = chamberLevelObjects.seedMatching[seed];
m_out_SP->push_back(*seed->parentBucket());
fillSeedInfo(seed, seedMatch.matchProb);
@@ -300,7 +307,7 @@ namespace MuonValR4 {
}
// also fill the reco not matched to any truth
for (auto & [ seed, assocInfo ] : chamberLevelObjects.seedMatching) {
- if (assocInfo.truthParticle) continue;
+ if (assocInfo.truthsegment) continue;
fillChamberInfo(chamber);
m_out_SP->push_back(*seed->parentBucket());
fillSeedInfo(seed, 0.);
diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.h b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.h
index 69109b96a1a3ad3cafa0a3f30f9df39ab4d9712b..103ff3490b1f13cf98fbeaa04306544f8ac9b326 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonHoughTransformTester.h
@@ -8,11 +8,13 @@
// Framework includes
#include "AthenaBaseComps/AthHistogramAlgorithm.h"
#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/ReadDecorHandle.h"
#include "StoreGate/ReadHandleKeyArray.h"
#include "StoreGate/ReadCondHandleKey.h"
// EDM includes
#include "xAODMuonSimHit/MuonSimHitContainer.h"
+#include "xAODMuon/MuonSegmentContainer.h"
#include <MuonPatternEvent/MuonPatternContainer.h>
@@ -60,7 +62,7 @@ namespace MuonValR4{
struct chamberLevelObjects {
struct SeedMatchQuantites {
/** @brief Best matched truth particle */
- HepMC::ConstGenParticlePtr truthParticle{};
+ const xAOD::MuonSegment* truthsegment{nullptr};
/** @brief Probability of which the segment is matched to it */
double matchProb{0.};
/** @brief Associated segment */
@@ -71,25 +73,25 @@ namespace MuonValR4{
/** @brief Collection of the truth particle trajectory */
struct TruthMatchQuantities{
- std::vector<const xAOD::MuonSimHit*> detectorHits{};
+ const xAOD::MuonSegment* truthsegment{nullptr};
std::vector<const MuonR4::SegmentSeed*> assocSeeds{};
};
- std::map<HepMC::ConstGenParticlePtr, TruthMatchQuantities> truthMatching{};
+ std::map<const xAOD::MuonSegment*, TruthMatchQuantities> truthMatching{};
};
- void matchSeedToTruth(const MuonR4::SegmentSeed* seed, chamberLevelObjects & objs ) const;
- std::pair<HepMC::ConstGenParticlePtr, double> matchSegmentToTruth(const MuonR4::Segment* seed, chamberLevelObjects & objs ) const;
- void matchSeedsToTruth(chamberLevelObjects & objs) const;
- void matchSegmentsToTruth(chamberLevelObjects & objs) const;
+ void matchSeedToTruth(const EventContext & ctx, const MuonR4::SegmentSeed* seed, chamberLevelObjects & objs ) const;
+ std::pair<const xAOD::MuonSegment*, double> matchSegmentToTruth(const MuonR4::Segment* seed, chamberLevelObjects & objs ) const;
+ void matchSeedsToTruth(const EventContext & ctx,chamberLevelObjects & objs) const;
void fillChamberInfo(const MuonGMR4::MuonChamber* chamber);
- void fillTruthInfo(const HepMC::ConstGenParticlePtr genParticlePtr, const std::vector<const xAOD::MuonSimHit*> & simHits, const ActsGeometryContext & gctx);
+ void fillTruthInfo(const EventContext & ctx, const xAOD::MuonSegment* truthSegment, const ActsGeometryContext & gctx);
void fillSeedInfo(const MuonR4::SegmentSeed* segmentSeed, double matchProb);
void fillSegmentInfo(const ActsGeometryContext& gctx,
const MuonR4::Segment* segment, double matchProb);
// MDT sim hits in xAOD format
SG::ReadHandleKeyArray<xAOD::MuonSimHitContainer> m_inSimHitKeys {this, "SimHitKeys",{}, "xAOD SimHit collections"};
+ SG::ReadHandleKey<xAOD::MuonSegmentContainer> m_truthSegmentKey {this, "TruthSegmentKey","TruthSegmentsR4", "truth segment container"};
SG::ReadHandleKey<MuonR4::SegmentSeedContainer> m_inHoughSegmentSeedKey{this, "SegmentSeedKey", "MuonHoughStationSegmentSeeds"};
SG::ReadHandleKey<MuonR4::SegmentContainer> m_inSegmentKey{this, "SegmentKey", "R4MuonSegments"};
@@ -97,6 +99,10 @@ namespace MuonValR4{
SG::ReadHandleKey<ActsGeometryContext> m_geoCtxKey{this, "AlignmentKey", "ActsAlignment", "cond handle key"};
+ using SimHitLinkVec = std::vector<ElementLink<xAOD::MuonSimHitContainer>>;
+
+ SG::ReadDecorHandleKey<xAOD::MuonSegmentContainer> m_rh_truthSegmentSimHitLink{this, "TruthSegmentSimHitLink", m_truthSegmentKey, "simHitLinks"};
+
ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc{this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};
const MuonGMR4::MuonDetectorManager* m_r4DetMgr{nullptr};
@@ -123,7 +129,7 @@ namespace MuonValR4{
MuonVal::ScalarBranch<unsigned int>& m_out_gen_nMDTHits{m_tree.newScalar<unsigned int>("genNMdtHits",0)};
MuonVal::ScalarBranch<unsigned int>& m_out_gen_nTGCHits{m_tree.newScalar<unsigned int>("genNTgcHits",0)};
MuonVal::ScalarBranch<unsigned int>& m_out_gen_nsTGCHits{m_tree.newScalar<unsigned int>("genNsTgcHits",0)};
- MuonVal::ScalarBranch<unsigned int>& m_out_gen_nMMits{m_tree.newScalar<unsigned int>("genNMmHits",0)};
+ MuonVal::ScalarBranch<unsigned int>& m_out_gen_nNswHits{m_tree.newScalar<unsigned int>("genNNswHits",0)};
MuonVal::ScalarBranch<bool>& m_out_hasMax {m_tree.newScalar<bool>("hasMax", false)};
MuonVal::ScalarBranch<bool>& m_out_max_hasPhiExtension {m_tree.newScalar<bool>("maxHasPhiExtension", false)};
MuonVal::ScalarBranch<float>& m_out_max_matchFraction {m_tree.newScalar<float>("maxMatchFraction", false)};