208 {
212 return StatusCode::FAILURE;
213 }
214 SG::ReadHandle<TgcL0SegmentContainer> segments{
m_segmentKey, ctx};
217 return StatusCode::FAILURE;
218 }
222 return StatusCode::FAILURE;
223 }
224 auto output = std::make_unique<TgcL0ValidationEvent>();
225 output->event.runNumber = ctx.eventID().run_number();
226 output->event.eventNumber = ctx.eventID().event_number();
227 output->event.lumiBlock = ctx.eventID().lumi_block();
228 output->event.bcid = ctx.eventID().bunch_crossing_id();
229
233 if (readFinalCandidates) {
234 SG::ReadHandle<xAOD::TGCCandDataContainer> finalCandidates{
236 if (!finalCandidates.
isValid()) {
238 return StatusCode::FAILURE;
239 }
240 finalCandidateContainer = finalCandidates.
cptr();
241 for (std::size_t inputIndex = 0U;
242 inputIndex < finalCandidateContainer->
size(); ++inputIndex) {
244 (*finalCandidateContainer)[inputIndex];
245 if (inputCandidate == nullptr) {
246 ATH_MSG_ERROR(
"Null final TGC candidate at index " << inputIndex);
247 return StatusCode::FAILURE;
248 }
249
250 const float eta = decodeEta(*inputCandidate);
251 const float phi = decodePhi(*inputCandidate);
252 int sourceIndex{-1};
253 int station{-1};
254 if (inputCandidate->
tcId() != 0U) {
255 sourceIndex =
256 sourceCandidateIndex(*inputCandidate, *candidates,
eta,
phi);
257 if (sourceIndex < 0) {
259 << (sourceIndex == -2 ? " has ambiguous"
260 : " has no")
261 << " pre-Inner source candidate");
262 return StatusCode::FAILURE;
263 }
264 station = pivotStation(
265 (*candidates)[static_cast<std::size_t>(sourceIndex)]
266 .positionStationMask);
267 if (station < 0) {
269 << " has no valid reference station");
270 return StatusCode::FAILURE;
271 }
272 }
273
274 output->finalCandidates.sourceCandidateIndex.emplace_back(sourceIndex);
275 output->finalCandidates.referenceStation.emplace_back(
276 station < 0
278 : static_cast<std::uint8_t>(station));
279 output->finalCandidates.subdetectorId.emplace_back(
281 output->finalCandidates.triggerSector.emplace_back(
283 output->finalCandidates.bcTag.emplace_back(inputCandidate->
bcTag());
284 output->finalCandidates.tcId.emplace_back(inputCandidate->
tcId());
285 output->finalCandidates.rawEta.emplace_back(inputCandidate->
eta());
286 output->finalCandidates.rawPhi.emplace_back(inputCandidate->
phi());
287 output->finalCandidates.eta.emplace_back(
eta);
288 output->finalCandidates.phi.emplace_back(
phi);
289 output->finalCandidates.ptCode.emplace_back(inputCandidate->
pt());
291 output->finalCandidates.threshold.emplace_back(
293 output->finalCandidates.charge.emplace_back(
295 output->finalCandidates.innerCoincidence.emplace_back(
297 output->finalCandidates.goodMagneticField.emplace_back(
299 output->finalCandidates.truthIndex.emplace_back(-1);
300 output->finalCandidates.truthMatchDeltaR.emplace_back(
302 }
303 }
304
306 SG::ReadHandle<xAOD::SectorLogicCandDataContainer> sectorLogicCandidates{
308 if (!sectorLogicCandidates.
isValid()) {
310 return StatusCode::FAILURE;
311 }
312 if (!finalCandidateContainer)[[
unlikely]]{
314 return StatusCode::FAILURE;
315 }
316 std::size_t sectorLogicIndex{0U};
317 for (std::size_t inputIndex = 0U;
318 inputIndex < finalCandidateContainer->
size();
319 ++inputIndex) {
321 (*finalCandidateContainer)[inputIndex];
322 if (inputCandidate->
tcId() == 0U)
continue;
323 if (sectorLogicIndex >= sectorLogicCandidates->size()) {
325 "Fewer Sector Logic candidates than non-empty TGC candidates");
326 return StatusCode::FAILURE;
327 }
329 (*sectorLogicCandidates)[sectorLogicIndex];
330 if (sectorLogicCandidate == nullptr) {
332 << sectorLogicIndex);
333 return StatusCode::FAILURE;
334 }
335 if (!matchesTgcFields(*inputCandidate, *sectorLogicCandidate)) {
337 << " does not preserve TGC candidate " << inputIndex);
338 return StatusCode::FAILURE;
339 }
340 if (!hasOnlyTgcFields(*sectorLogicCandidate)) {
342 << " contains non-TGC payload");
343 return StatusCode::FAILURE;
344 }
345 if (sectorLogicCandidate->
boardID() != 0U ||
346 sectorLogicCandidate->
fiberID() != 0U ||
348 sectorLogicCandidate->
veto() != 0U) {
350 << " has unexpected placeholder metadata");
351 return StatusCode::FAILURE;
352 }
353 output->sectorLogic.inputCandidateIndex.emplace_back(
354 static_cast<std::uint32_t>(inputIndex));
355 output->sectorLogic.candWord.emplace_back(
357 output->sectorLogic.candExtraWord.emplace_back(
359 output->sectorLogic.boardId.emplace_back(sectorLogicCandidate->
boardID());
360 output->sectorLogic.fiberId.emplace_back(sectorLogicCandidate->
fiberID());
361 output->sectorLogic.bcidOffset.emplace_back(
363 output->sectorLogic.veto.emplace_back(sectorLogicCandidate->
veto());
364 ++sectorLogicIndex;
365 }
366 if (sectorLogicIndex != sectorLogicCandidates->size()) {
368 "More Sector Logic candidates than non-empty TGC candidates");
369 return StatusCode::FAILURE;
370 }
371 }
372
373 const auto processParticle = [&](
const auto&
particle) -> StatusCode {
374 if (!particle || std::abs(
particle->pdg_id()) != 13 ||
377 return StatusCode::SUCCESS;
378 }
379
386 return StatusCode::SUCCESS;
387 }
388
390 const double theta = 2.0 * std::atan(std::exp(-
eta));
391 const double momentumMagnitude =
pt * std::cosh(
eta);
392 if (momentumMagnitude <= 0.0 || !std::isfinite(
theta)) {
393 return StatusCode::SUCCESS;
394 }
395
397 const Trk::PerigeeSurface perigeeSurface{
Amg::Vector3D{0.0, 0.0, 0.0}};
399 charge / momentumMagnitude, perigeeSurface};
400
401 StationPositions positions{};
402 std::uint8_t stationMask{0U};
403 const std::vector<double>& stationAbsZ =
m_stationAbsZ.value();
404 for (std::size_t station = 0U; station < positions.size(); ++station) {
405 const double z =
eta >= 0.0 ? stationAbsZ[station] : -stationAbsZ[station];
410 const Trk::BoundaryCheck boundaryCheck{true};
417 continue;
418 }
419 positions[station] = StationPosition{
420 true,
static_cast<float>(
position.eta()),
422 stationMask |= stationBit(station);
423 }
424
425 output->truth.pdgId.emplace_back(pdgId);
427 output->truth.pt.emplace_back(
static_cast<float>(pt));
428 output->truth.eta.emplace_back(
static_cast<float>(
eta));
429 output->truth.phi.emplace_back(
static_cast<float>(
phi));
430 output->truth.charge.emplace_back(
static_cast<float>(
charge));
431 output->truth.extrapolatedStationMask.emplace_back(stationMask);
432 output->truth.m1Eta.emplace_back(positions[0].
eta);
433 output->truth.m1Phi.emplace_back(positions[0].
phi);
434 output->truth.m2Eta.emplace_back(positions[1].
eta);
435 output->truth.m2Phi.emplace_back(positions[1].
phi);
436 output->truth.m3Eta.emplace_back(positions[2].
eta);
437 output->truth.m3Phi.emplace_back(positions[2].
phi);
438 output->truth.matched.emplace_back(0U);
439 output->truth.matchedCandidateIndex.emplace_back(-1);
441 output->truth.unmatchedReason.emplace_back(
442 stationMask == 0x7U
443 ? static_cast<std::uint8_t>(
445 : static_cast<std::uint8_t>(
447 output->truth.finalCandidateMatched.emplace_back(0U);
448 output->truth.matchedFinalCandidateIndex.emplace_back(-1);
449 output->truth.finalCandidateMatchDeltaR.emplace_back(
451 output->truth.finalCandidateUnmatchedReason.emplace_back(
452 stationMask == 0x7U
453 ? static_cast<std::uint8_t>(
455 : static_cast<std::uint8_t>(
457 output->truth.wireSegmentMatched.emplace_back(0U);
458 output->truth.stripSegmentMatched.emplace_back(0U);
459 output->truth.matchedWireSegmentIndex.emplace_back(-1);
460 output->truth.matchedStripSegmentIndex.emplace_back(-1);
461 output->truth.wireSegmentMatchResidual.emplace_back(
463 output->truth.stripSegmentMatchResidual.emplace_back(
465 return StatusCode::SUCCESS;
466 };
467
469 if (event == nullptr) continue;
470#if __has_include("HepMC3/GenEvent.h")
471 for (
const auto& particle :
event->particles()) {
473 }
474#else
475 for (
auto particle =
event->particles_begin();
476 particle !=
event->particles_end(); ++particle) {
478 }
479#endif
480 }
481
482 for (const TgcL0Candidate& candidate : *candidates) {
483 output->candidates.subdetectorId.emplace_back(candidate.subdetectorId);
484 output->candidates.triggerSector.emplace_back(candidate.sectorId);
485 output->candidates.readoutSector.emplace_back(candidate.readoutSector);
486 output->candidates.bcTag.emplace_back(candidate.bcTag);
487 output->candidates.stationMask.emplace_back(candidate.positionStationMask);
488 output->candidates.wireStationMask.emplace_back(candidate.wireStationMask);
489 output->candidates.stripStationMask.emplace_back(candidate.stripStationMask);
490 output->candidates.eta.emplace_back(candidate.eta);
491 output->candidates.phi.emplace_back(candidate.phi);
492 output->candidates.deltaTheta.emplace_back(candidate.deltaTheta);
493 output->candidates.deltaPhi.emplace_back(candidate.deltaPhi);
494 output->candidates.pt.emplace_back(candidate.preInnerCoincidencePt);
495 output->candidates.threshold.emplace_back(
496 candidate.preInnerCoincidenceThreshold);
497 output->candidates.charge.emplace_back(candidate.charge);
498 output->candidates.goodMagneticField.emplace_back(
499 candidate.goodMagneticField ? 1U : 0U);
500 output->candidates.truthIndex.emplace_back(-1);
501 }
502
503 std::vector<CandidateMatch> candidateMatches;
504 std::vector<bool> hasCandidateInWindow(
output->truth.pt.size(),
false);
505 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
506 if (
output->truth.extrapolatedStationMask[truth] != 0x7U)
continue;
507 const StationPositions positions = truthPositions(*output, truth);
508 for (std::size_t candidate = 0U; candidate <
candidates->size();
509 ++candidate) {
510 const TgcL0Candidate& inputCandidate = (*candidates)[candidate];
513 continue;
514 }
515 if (
output->truth.eta[truth] * inputCandidate.eta < 0.F)
continue;
516 const float residual = meanDeltaR(inputCandidate, positions);
517 if (!std::isfinite(residual) || residual >
m_maxMeanDeltaR.value())
continue;
518 hasCandidateInWindow[truth] = true;
519 candidateMatches.push_back({
residual, truth, candidate});
520 }
521 }
523 [](const CandidateMatch& lhs, const CandidateMatch& rhs) {
524 if (lhs.residual != rhs.residual) {
525 return lhs.residual < rhs.residual;
526 }
527 if (lhs.truth != rhs.truth)
return lhs.truth < rhs.truth;
528 return lhs.candidate < rhs.candidate;
529 });
530 for (
const CandidateMatch&
match : candidateMatches) {
531 if (
output->truth.matchedCandidateIndex[
match.truth] >= 0 ||
532 output->candidates.truthIndex[
match.candidate] >= 0) {
533 continue;
534 }
537 static_cast<int>(
match.candidate);
542 static_cast<int>(
match.truth);
543 }
544 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
545 if (
output->truth.matched[truth] != 0U ||
546 output->truth.extrapolatedStationMask[truth] != 0x7U) {
547 continue;
548 }
549 output->truth.unmatchedReason[truth] =
static_cast<std::uint8_t
>(
550 hasCandidateInWindow[truth]
553 }
554
556 std::vector<FinalCandidateMatch> finalCandidateMatches;
557 std::vector<bool> hasFinalCandidateInWindow(
output->truth.pt.size(),
558 false);
559 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
560 if (
output->truth.extrapolatedStationMask[truth] != 0x7U)
continue;
561 const StationPositions positions = truthPositions(*output, truth);
562 for (std::size_t candidate = 0U;
563 candidate <
output->finalCandidates.tcId.size(); ++candidate) {
564 if (
output->finalCandidates.tcId[candidate] == 0U)
continue;
566 (
output->finalCandidates.bcTag[candidate] &
568 continue;
569 }
570 if (
output->truth.eta[truth] *
571 output->finalCandidates.eta[candidate] <
572 0.F) {
573 continue;
574 }
575 const std::size_t station =
576 output->finalCandidates.referenceStation[candidate];
577 if (station >= positions.size() || !positions[station].valid) continue;
579 positions[station].eta;
580 const float deltaPhi =
static_cast<float>(
582 positions[station].phi));
585 continue;
586 }
587 hasFinalCandidateInWindow[truth] = true;
588 finalCandidateMatches.push_back({
residual, truth, candidate});
589 }
590 }
592 finalCandidateMatches.begin(), finalCandidateMatches.end(),
593 [](const FinalCandidateMatch& lhs, const FinalCandidateMatch& rhs) {
594 if (lhs.residual != rhs.residual) {
595 return lhs.residual < rhs.residual;
596 }
597 if (lhs.truth != rhs.truth)
return lhs.truth < rhs.truth;
598 return lhs.candidate < rhs.candidate;
599 });
600 for (
const FinalCandidateMatch&
match : finalCandidateMatches) {
601 if (
output->truth.matchedFinalCandidateIndex[
match.truth] >= 0 ||
602 output->finalCandidates.truthIndex[
match.candidate] >= 0) {
603 continue;
604 }
605 output->truth.finalCandidateMatched[
match.truth] = 1U;
606 output->truth.matchedFinalCandidateIndex[
match.truth] =
607 static_cast<int>(
match.candidate);
609 output->truth.finalCandidateUnmatchedReason[
match.truth] =
611 output->finalCandidates.truthIndex[
match.candidate] =
612 static_cast<int>(
match.truth);
613 output->finalCandidates.truthMatchDeltaR[
match.candidate] =
615 }
616 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
617 if (
output->truth.finalCandidateMatched[truth] != 0U ||
618 output->truth.extrapolatedStationMask[truth] != 0x7U) {
619 continue;
620 }
621 output->truth.finalCandidateUnmatchedReason[truth] =
622 static_cast<std::uint8_t>(
623 hasFinalCandidateInWindow[truth]
626 }
627 }
628
629 for (const TgcL0Segment& segment : *segments) {
630 output->segments.subdetectorId.emplace_back(segment.subdetectorId);
631 output->segments.triggerSector.emplace_back(segment.triggerSector);
632 output->segments.bcTag.emplace_back(segment.bcTag);
633 output->segments.projection.emplace_back(
634 static_cast<std::uint8_t>(segment.projection));
635 output->segments.stationMask.emplace_back(segment.stationMask);
636 output->segments.summedQuality.emplace_back(segment.summedQuality);
637 output->segments.nStations.emplace_back(segment.nStations);
638 output->segments.eta.emplace_back(segment.eta);
639 output->segments.phi.emplace_back(segment.phi);
640 output->segments.residual.emplace_back(segment.residual);
641 output->segments.outputResidual.emplace_back(segment.outputResidual);
642 output->segments.consistency.emplace_back(segment.consistency);
643 output->segments.pivotChannel.emplace_back(segment.pivotChannel);
644 output->segments.truthIndex.emplace_back(-1);
645 output->segments.truthMatchResidual.emplace_back(
647 }
648
650 const float maximumResidual) {
651 const auto projectionValue = static_cast<std::uint8_t>(projection);
652 std::vector<SegmentMatch> possibleMatches;
653 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
654 for (std::size_t segment = 0U;
655 segment <
output->segments.projection.size(); ++segment) {
656 if (
output->segments.projection[segment] != projectionValue)
continue;
658 (
output->segments.bcTag[segment] &
660 continue;
661 }
662 if (
output->truth.eta[truth] *
output->segments.eta[segment] < 0.F) {
663 continue;
664 }
665 const int station = pivotStation(
output->segments.stationMask[segment]);
666 if (station < 0 ||
667 (
output->truth.extrapolatedStationMask[truth] &
668 stationBit(static_cast<std::size_t>(station))) == 0U) {
669 continue;
670 }
671 const std::array<float, 3>
truthEta{
673 output->truth.m3Eta[truth]};
674 const std::array<float, 3>
truthPhi{
676 output->truth.m3Phi[truth]};
679 ? std::abs(
output->segments.eta[segment] - truthEta[station])
682 if (std::isfinite(residual) && residual <= maximumResidual) {
683 possibleMatches.push_back({
residual, truth, segment});
684 }
685 }
686 }
688 [](const SegmentMatch& lhs, const SegmentMatch& rhs) {
689 if (lhs.residual != rhs.residual) {
690 return lhs.residual < rhs.residual;
691 }
692 if (lhs.truth != rhs.truth)
return lhs.truth < rhs.truth;
693 return lhs.segment < rhs.segment;
694 });
695 for (
const SegmentMatch&
match : possibleMatches) {
696 int& truthSegment =
698 ?
output->truth.matchedWireSegmentIndex[
match.truth]
699 :
output->truth.matchedStripSegmentIndex[
match.truth];
700 if (truthSegment >= 0 ||
702 continue;
703 }
704 truthSegment =
static_cast<int>(
match.segment);
706 static_cast<int>(
match.truth);
709 output->truth.wireSegmentMatched[
match.truth] = 1U;
711 } else {
712 output->truth.stripSegmentMatched[
match.truth] = 1U;
714 }
715 }
716 };
717
722
725 ATH_MSG_ERROR(
"Refusing to record inconsistent validation data: "
727 return StatusCode::FAILURE;
728 }
729
730 SG::WriteHandle<TgcL0ValidationEvent> outputHandle{
m_outputKey, ctx};
732 return StatusCode::SUCCESS;
733}
Scalar eta() const
pseudorapidity method
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x,...)
double charge(const T &p)
size_type size() const noexcept
Returns the number of elements in the collection.
Gaudi::Property< float > m_maxMeanDeltaR
Gaudi::Property< int > m_requiredTruthStatus
SG::WriteHandleKey< TgcL0ValidationEvent > m_outputKey
Gaudi::Property< std::uint16_t > m_requiredBcTagMask
Gaudi::Property< double > m_minPt
Gaudi::Property< std::vector< double > > m_stationAbsZ
Gaudi::Property< int > m_maxAbsBarcode
Gaudi::Property< double > m_validationMinR
SG::ReadHandleKey< TgcL0CandidateContainer > m_candidateKey
SG::ReadHandleKey< xAOD::TGCCandDataContainer > m_finalCandidateKey
Gaudi::Property< bool > m_validateSectorLogic
Gaudi::Property< double > m_validationMaxR
SG::ReadHandleKey< TgcL0SegmentContainer > m_segmentKey
Gaudi::Property< float > m_maxStripSegmentDeltaPhi
Gaudi::Property< double > m_minAbsEta
Gaudi::Property< bool > m_validateFinalCandidates
Gaudi::Property< double > m_maxAbsEta
Gaudi::Property< float > m_maxWireSegmentDeltaEta
Gaudi::Property< float > m_maxFinalCandidateDeltaR
Gaudi::Property< double > m_validationPlaneToleranceZ
ToolHandle< Trk::IExtrapolator > m_extrapolator
SG::ReadHandleKey< xAOD::SectorLogicCandDataContainer > m_sectorLogicKey
SG::ReadHandleKey< McEventCollection > m_truthEventKey
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
uint16_t sectorId() const
Retrieve the sector id.
uint16_t eta() const
Retrieve the eta.
uint8_t candCharge() const
Retrieve the candidate charge.
uint8_t threshold() const
Retrieve the threshold.
uint16_t bcTag() const
Retrieve the bunch crossing tag.
uint8_t pt() const
Retrieve the encoded candidate pT.
uint16_t phi() const
Retrieve the phi.
uint16_t subdetectorId() const
Retrieve the sub detector id.
uint16_t fiberID() const
Retrieve the fiber ID.
unsigned short veto() const
Retrieve the veto flag.
uint32_t candWord() const
Retrieve the candidate word (first 32-bits).
int BCIDOffset() const
Retrieve the bunch crossing identifier offset.
uint16_t boardID() const
Retrieve the board ID.
uint32_t candExtraWord() const
Retrieve the candidate extra word (second 32-bits).
bool goodMagneticField() const
Check whether the candidate is in a good magnetic-field region.
bool hasInnerCoincidence() const
Check whether the candidate passed Inner Coincidence.
float ptValueGeV() const
Retrieve the inherited encoded pT value in GeV.
uint8_t tcId() const
Retrieve the trigger-candidate identifier.
bool match(std::string s1, std::string s2)
match the individual directories of two strings
Eigen::Affine3d Transform3D
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Eigen::Matrix< double, 3, 1 > Vector3D
HepMC3::GenEvent GenEvent
TgcL0ValidationProjection
Projection encoding used in the validation segment block.
TgcL0ValidationCheckResult checkTgcL0ValidationEvent(const TgcL0ValidationEvent &event)
Check vector sizes and reciprocal cross-block indices.
static constexpr float TgcL0ValidationInvalidValue
Common sentinel for unavailable floating-point validation data.
static constexpr std::uint8_t TgcL0ValidationNoUnmatchedReason
Sentinel indicating that a matched truth object has no failure code.
double deltaEta(const I4Momentum &p1, const I4Momentum &p2)
Computes efficiently .
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
void stable_sort(DataModel_detail::iterator< DVL > beg, DataModel_detail::iterator< DVL > end)
Specialization of stable_sort for DataVector/List.
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
TGCCandDataContainer_v1 TGCCandDataContainer
SectorLogicCandData_v1 SectorLogicCandData
TGCCandData_v1 TGCCandData