6 #include "GaudiKernel/SystemOfUnits.h"
14 #include "Acts/Utilities/Enumerate.hpp"
15 #include "GaudiKernel/PhysicalConstants.h"
29 using simHitSet = std::unordered_set<const xAOD::MuonSimHit*>;
47 template <
class SpType>
60 m_tree.addBranch(std::make_unique<EventInfoBranch>(m_tree, infoOpts));
63 ATH_CHECK(m_truthSegmentKey.initialize(!m_truthSegmentKey.empty()));
66 ATH_CHECK(m_inHoughSegmentSeedKeys.initialize());
67 ATH_CHECK(m_spKey.initialize(m_writeSpacePoints));
68 if (m_writeSpacePoints) {
69 m_spTester = std::make_unique<SpacePointTesterModule>(m_tree, m_spKey.key(), msgLevel());
70 m_tree.addBranch(m_spTester);
72 m_tree.disableBranch(m_spMatchedToPattern.name());
73 m_tree.disableBranch(m_spMatchedToSegment.name());
79 ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
81 return StatusCode::SUCCESS;
89 using namespace SegmentFit;
92 const std::vector<int> truthSigns = SeedingAux::strawSigns(truePos, trueDir, recoSeg.
measurements());
93 const std::vector<int> recoSigns = SeedingAux::strawSigns(recoPos, recoDir, recoSeg.
measurements());
94 for (
unsigned int s = 0 ;
s < truthSigns.size(); ++
s) {
95 same += (truthSigns[
s] != 0) && truthSigns[
s] == recoSigns[
s];
99 std::vector<ObjectMatching>
104 std::vector<ObjectMatching> allAssociations{};
105 std::unordered_set<const SegmentSeed*> usedSeeds{};
106 std::unordered_set<const Segment*> usedSegs{};
112 std::vector<simHitSet> truthHitsVec{}, seedSimHitVec{}, segmentSimHitVec{};
116 for (
const Segment* segment: *segmentContainer){
124 ObjectMatching & matchedWithTruth = allAssociations.emplace_back();
126 matchedWithTruth.
chamber = m_detMgr->getSectorEnvelope((*truthHits.begin())->identify());
128 std::vector<std::pair<const SegmentSeed*, unsigned>> matchedSeeds{};
135 if (seed->msSector() != matchedWithTruth.
chamber) {
138 const simHitSet& seedHits{seedSimHitVec[seedIdx]};
139 unsigned int matchedHits =
countMatched(truthHits, seedHits);
143 matchedSeeds.emplace_back(std::make_pair(seed, matchedHits));
146 std::vector<std::pair<const Segment*, unsigned>> matchedSegs{};
148 for (
const Segment* segment : *segmentContainer) {
150 if (segment->msSector() != matchedWithTruth.
chamber) {
153 const simHitSet& segmentHits{segmentSimHitVec[segmentIdx]};
154 unsigned int matchedHits =
countMatched(truthHits, segmentHits);
158 matchedSegs.emplace_back(std::make_pair(segment,matchedHits));
164 std::ranges::sort(matchedSegs,
165 [
this, &truth, &gctx](
const std::pair<const Segment*, unsigned>& segA,
166 const std::pair<const Segment*, unsigned>& segB){
167 if (segA.second != segB.second)
return segA.second > segB.second;
168 return countOnSameSide(gctx, *truth, *segA.first) > countOnSameSide(gctx, *truth, *segB.first);
171 std::ranges::sort(matchedSeeds, [](
const std::pair<const SegmentSeed*, unsigned>& seedA,
172 const std::pair<const SegmentSeed*, unsigned>& seedB) {
173 return seedA.second > seedB.second;
180 for (
const auto& [
matched, nMatchedHits] : matchedSegs) {
184 usedSeeds.insert(
matched->parent());
189 for (
const auto& [seed , nHits] : matchedSeeds) {
194 usedSeeds.insert(seed);
204 for (
const Segment* seg: *segmentContainer) {
206 if (usedSegs.count(seg)) {
210 match.chamber = seg->msSector();
211 match.matchedSegments = {seg};
212 match.matchedSeeds = {seg->parent()};
214 usedSeeds.insert(seg->parent());
215 match.matchedSeedFoundSegment.push_back(1);
219 if (usedSeeds.count(seed)) {
223 match.chamber = seed->msSector();
224 match.matchedSeeds = {seed};
225 match.matchedSeedFoundSegment.push_back(0);
227 return allAssociations;
232 return StatusCode::SUCCESS;
236 const EventContext & ctx = Gaudi::Hive::currentContext();
248 segmentSeeds.insert(segmentSeeds.end(),readSegmentSeeds->begin(), readSegmentSeeds->end());
253 segments.insert(segments.end(),readSegments->begin(), readSegments->end());
261 ATH_MSG_DEBUG(
"Succesfully retrieved input collections. Seeds: "<<segmentSeeds.size()
262 <<
", segments: "<<segments.size() <<
", truth segments: "<<(readTruthSegments? readTruthSegments->size() : -1)<<
".");
263 std::vector<ObjectMatching>
objects = matchWithTruth(gctx, readTruthSegments, segmentSeeds.asDataVector(),
264 segments.asDataVector());
266 fillChamberInfo(
obj.chamber);
268 fillSegmentInfo(gctx,
obj);
269 if(m_isMC) fillTruthInfo(gctx,
obj.truthSegment);
272 return StatusCode::SUCCESS;
275 m_out_chamberIndex = Acts::toUnderlying(msSector->
chamberIndex());
276 m_out_stationSide = msSector->
side();
281 if (!segment)
return;
282 m_out_hasTruth =
true;
288 m_out_gen_Eta = segDir.eta();
289 m_out_gen_Phi = segDir.phi();
290 m_out_gen_Pt = acc_pt(*segment);
291 m_out_gen_Q = acc_charge(*segment);
303 m_out_gen_y0 = chamberPos.y();
304 m_out_gen_x0 = chamberPos.x();
305 m_out_gen_time = segment->
t0();
315 const Amg::Vector3D chamberPos = localToChamber * xAOD::toEigen(hit->localPosition());
316 minYhit =
std::min(chamberPos.y(), minYhit);
317 maxYhit =
std::max(chamberPos.y(), maxYhit);
319 m_out_gen_minYhit = minYhit;
320 m_out_gen_maxYhit = maxYhit;
322 ATH_MSG_DEBUG(
"A true max on chamber index "<<m_out_chamberIndex.getVariable()<<
" side "<<m_out_stationSide.getVariable()<<
" phi "<<m_out_stationPhi.getVariable()<<
" with "
323 <<m_out_gen_nMDTHits.getVariable()<<
" MDT and "<<m_out_gen_nRPCHits.getVariable()+m_out_gen_nTGCHits.getVariable()<<
" trigger hits is at "
324 <<m_out_gen_tantheta.getVariable()<<
" and "<<m_out_gen_y0.getVariable());
329 m_out_nSpacePoints = bucket.size();
330 m_out_nPrecSpacePoints = std::ranges::count_if(bucket, [](
const SpacePointBucket::value_type& sp){
333 m_out_nPhiSpacePoints = std::ranges::count_if(bucket, [](
const SpacePointBucket::value_type& sp){
334 return sp->measuresPhi();
336 if (!m_visionTool.isEnabled()){
339 m_out_nTrueSpacePoints = std::ranges::count_if(bucket,[
this](
const SpacePointBucket::value_type& sp){
340 return m_visionTool->isLabeled(*sp);
342 m_out_nTruePrecSpacePoints = std::ranges::count_if(bucket,[
this](
const SpacePointBucket::value_type& sp){
343 return isPrecHit(*sp) && m_visionTool->isLabeled(*sp);
345 m_out_nTruePhiSpacePoints = std::ranges::count_if(bucket,[
this](
const SpacePointBucket::value_type& sp){
346 return sp->measuresPhi() && m_visionTool->isLabeled(*sp);
351 m_out_seed_n =
obj.matchedSeeds.size();
352 for (
const auto [iseed, seed] : Acts::enumerate(
obj.matchedSeeds)){
354 fillBucketInfo(*seed->parentBucket());
356 double minYhit = m_out_bucketEnd.getVariable();
357 double maxYhit = m_out_bucketStart.getVariable();
358 for (
const SpacePoint* hit : seed->getHitsInMax()){
359 minYhit =
std::min(hit->localPosition().y(),minYhit);
360 maxYhit =
std::max(hit->localPosition().y(),maxYhit);
362 m_out_seed_minYhit.push_back(minYhit);
363 m_out_seed_maxYhit.push_back(maxYhit);
365 m_out_seed_hasPhiExtension.push_back(seed->hasPhiExtension());
366 m_out_seed_nMatchedHits.push_back(
countMatched(
obj.truthSegment, seed));
367 m_out_seed_y0.push_back(seed->interceptY());
368 m_out_seed_tantheta.push_back(seed->tanBeta());
369 if (seed->hasPhiExtension()){
370 m_out_seed_x0.push_back(seed->interceptX());
371 m_out_seed_tanphi.push_back(seed->tanAlpha());
373 m_out_seed_x0.push_back(-999);
374 m_out_seed_tanphi.push_back(-999);
377 m_out_seed_nHits.push_back(seed->getHitsInMax().size());
378 unsigned nMdtSeed{0}, nRpcSeed{0}, nTgcSeed{0}, nMmSeed{0}, nsTgcSeed{0};
379 unsigned nPrecHits{0}, nEtaHits{0}, nPhiHits{0}, nTrueHits{0}, nTruePrecHits{0}, nTrueEtaHits{0}, nTruePhiHits{0};
380 std::vector<unsigned char> treeIdxs{};
382 for (
const HoughHitType & houghSP: seed->getHitsInMax()){
383 if (m_writeSpacePoints){
384 unsigned treeIdx = m_spTester->push_back(*houghSP);
385 treeIdxs.push_back(treeIdx);
388 nPhiHits += houghSP->measuresPhi();
389 nEtaHits += houghSP->measuresEta();
391 if (m_visionTool.isEnabled()) {
392 nTrueHits += m_visionTool->isLabeled(*houghSP);
393 nTruePrecHits += m_visionTool->isLabeled(*houghSP) &&
isPrecHit(*houghSP);
394 nTruePhiHits += m_visionTool->isLabeled(*houghSP) && houghSP->measuresPhi();
395 nTrueEtaHits += m_visionTool->isLabeled(*houghSP) && houghSP->measuresEta();
397 switch (houghSP->type()) {
402 nRpcSeed+=houghSP->measuresEta();
403 nRpcSeed+=houghSP->measuresPhi();
406 nTgcSeed+=houghSP->measuresEta();
407 nTgcSeed+=houghSP->measuresPhi();
410 nsTgcSeed += houghSP->measuresEta();
411 nsTgcSeed += houghSP->measuresPhi();
417 ATH_MSG_WARNING(
"Technology "<<houghSP->identify() <<
" not yet implemented");
420 m_out_seed_nMdt.push_back(nMdtSeed);
421 m_out_seed_nRpc.push_back(nRpcSeed);
422 m_out_seed_nTgc.push_back(nTgcSeed);
423 m_out_seed_nsTgc.push_back(nsTgcSeed);
424 m_out_seed_nMm.push_back(nMmSeed);
426 m_out_seed_nPrecHits.push_back(nPrecHits);
427 m_out_seed_nEtaHits.push_back(nEtaHits);
428 m_out_seed_nPhiHits.push_back(nPhiHits);
430 m_out_seed_nTrueHits.push_back(nTrueHits);
431 m_out_seed_nTruePrecHits.push_back(nTruePrecHits);
432 m_out_seed_nTrueEtaHits.push_back(nTrueEtaHits);
433 m_out_seed_nTruePhiHits.push_back(nTruePhiHits);
435 m_out_seed_ledToSegment.push_back(
obj.matchedSeedFoundSegment.at(iseed));
436 if (m_writeSpacePoints) {
437 m_spMatchedToPattern[iseed] = std::move(treeIdxs);
445 using namespace SegmentFit;
447 m_out_segment_n =
obj.matchedSegments.size();
448 for (
auto & segment :
obj.matchedSegments){
449 m_out_segment_hasPhi.push_back(std::ranges::find_if(segment->measurements(), [](
const auto& meas){ return meas->measuresPhi();})
450 !=segment->measurements().end());
451 m_out_segment_fitIter.push_back(segment->nFitIterations());
452 m_out_segment_truthMatchedHits.push_back(
countMatched(
obj.truthSegment, segment));
453 m_out_segment_chi2.push_back(segment->chi2());
454 m_out_segment_nDoF.push_back(segment->nDoF());
455 m_out_segment_hasTimeFit.push_back(segment->hasTimeFit());
457 m_out_segment_err_x0.push_back(segment->covariance()(Acts::toUnderlying(ParamDefs::x0), Acts::toUnderlying(ParamDefs::x0)));
458 m_out_segment_err_y0.push_back(segment->covariance()(Acts::toUnderlying(ParamDefs::y0), Acts::toUnderlying(ParamDefs::y0)));
461 m_out_segment_err_time.push_back(segment->covariance()(Acts::toUnderlying(
ParamDefs::t0), Acts::toUnderlying(
ParamDefs::t0)));
465 m_out_segment_y0.push_back(locPos.y());
466 m_out_segment_x0.push_back(locPos.x());
467 m_out_segment_time.push_back(segment->segementT0() + segment->position().mag() * c_inv);
469 unsigned nMdtHits{0}, nRpcEtaHits{0}, nRpcPhiHits{0}, nTgcEtaHits{0}, nTgcPhiHits{0},
470 nMmEtaHits{0}, nMmStereoHits{0}, nStgcStripHits{0},nStgcWireHits{0}, nStgcPadHits{0};
471 unsigned nTrueHits{0}, nTruePrecHits{0}, nTrueEtaHits{0}, nTruePhiHits{0};
476 std::vector<unsigned char>
matched;
477 for (
const auto & meas : segment->measurements()){
480 minYhit =
std::min(meas->localPosition().y(),minYhit);
481 maxYhit =
std::max(meas->localPosition().y(),maxYhit);
482 if (m_writeSpacePoints) {
483 unsigned treeIdx = m_spTester->push_back(*meas->spacePoint());
484 if (treeIdx >=
matched.size()){
489 if (m_visionTool.isEnabled()) {
490 nTrueHits += m_visionTool->isLabeled(*meas->spacePoint());
491 nTruePrecHits +=
isPrecHit(*meas) && m_visionTool->isLabeled(*meas->spacePoint());
492 nTrueEtaHits += meas->measuresEta() && m_visionTool->isLabeled(*meas->spacePoint());
493 nTruePhiHits += meas->measuresPhi() && m_visionTool->isLabeled(*meas->spacePoint());
495 switch (meas->type()) {
500 nRpcEtaHits += meas->measuresEta();
501 nRpcPhiHits += meas->measuresPhi();
504 nTgcEtaHits += meas->measuresEta();
505 nTgcPhiHits += meas->measuresPhi();
508 const MmIdHelper& idHelper{m_idHelperSvc->mmIdHelper()};
509 nMmEtaHits += !idHelper.
isStereo(meas->spacePoint()->identify());
510 nMmStereoHits += !idHelper.isStereo(meas->spacePoint()->identify());
513 const auto* prd =
static_cast<const xAOD::sTgcMeasurement*
>(meas->spacePoint()->primaryMeasurement());
514 nStgcStripHits += prd->
channelType() == sTgcIdHelper::sTgcChannelTypes::Strip;
515 nStgcWireHits += prd->channelType() == sTgcIdHelper::sTgcChannelTypes::Wire;
516 nStgcPadHits += prd->channelType() == sTgcIdHelper::sTgcChannelTypes::Pad;
519 nStgcWireHits += prd->
channelType() == sTgcIdHelper::sTgcChannelTypes::Wire;
520 nStgcPadHits += prd->channelType() == sTgcIdHelper::sTgcChannelTypes::Pad;
527 m_out_segment_nMdtHits.push_back(nMdtHits);
528 m_out_segment_nRpcEtaHits.push_back(nRpcEtaHits);
529 m_out_segment_nRpcPhiHits.push_back(nRpcPhiHits);
530 m_out_segment_nTgcEtaHits.push_back(nTgcEtaHits);
531 m_out_segment_nTgcPhiHits.push_back(nTgcPhiHits);
533 m_out_segment_nMmEtaHits.push_back(nMmEtaHits);
534 m_out_segment_nMmStereoHits.push_back(nMmStereoHits);
535 m_out_segment_nsTgcStripHits.push_back(nStgcStripHits);
536 m_out_segment_nsTgcWireHits.push_back(nStgcWireHits);
537 m_out_segment_nsTgcPadpHits.push_back(nStgcPadHits);
540 m_out_segment_nTrueHits.push_back(nTrueHits);
541 m_out_segment_nTruePrecHits.push_back(nTruePrecHits);
542 m_out_segment_nTruePhiHits.push_back(nTruePhiHits);
543 m_out_segment_nTrueEtaHits.push_back(nTrueEtaHits);
545 m_out_segment_minYhit.push_back(minYhit);
546 m_out_segment_maxYhit.push_back(maxYhit);
547 if (m_writeSpacePoints) {
548 m_spMatchedToSegment.push_back(std::move(
matched));