145 {
149 return StatusCode::FAILURE;
150 }
151 SG::ReadHandle<TgcL0SegmentContainer> segments{
m_segmentKey, ctx};
154 return StatusCode::FAILURE;
155 }
159 return StatusCode::FAILURE;
160 }
161
162 auto output = std::make_unique<TgcL0ValidationEvent>();
163 output->event.runNumber = ctx.eventID().run_number();
164 output->event.eventNumber = ctx.eventID().event_number();
165 output->event.lumiBlock = ctx.eventID().lumi_block();
166 output->event.bcid = ctx.eventID().bunch_crossing_id();
167
168 const auto processParticle = [&](
const auto&
particle) -> StatusCode {
169 if (!particle || std::abs(
particle->pdg_id()) != 13 ||
172 return StatusCode::SUCCESS;
173 }
174
181 return StatusCode::SUCCESS;
182 }
183
185 const double theta = 2.0 * std::atan(std::exp(-
eta));
186 const double momentumMagnitude =
pt * std::cosh(
eta);
187 if (momentumMagnitude <= 0.0 || !std::isfinite(
theta)) {
188 return StatusCode::SUCCESS;
189 }
190
192 const Trk::PerigeeSurface perigeeSurface{
Amg::Vector3D{0.0, 0.0, 0.0}};
194 charge / momentumMagnitude, perigeeSurface};
195
196 StationPositions positions{};
197 std::uint8_t stationMask{0U};
198 const std::vector<double>& stationAbsZ =
m_stationAbsZ.value();
199 for (std::size_t station = 0U; station < positions.size(); ++station) {
200 const double z =
eta >= 0.0 ? stationAbsZ[station] : -stationAbsZ[station];
205 const Trk::BoundaryCheck boundaryCheck{true};
212 continue;
213 }
214 positions[station] = StationPosition{
215 true,
static_cast<float>(
position.eta()),
217 stationMask |= stationBit(station);
218 }
219
220 output->truth.pdgId.emplace_back(pdgId);
222 output->truth.pt.emplace_back(
static_cast<float>(pt));
223 output->truth.eta.emplace_back(
static_cast<float>(
eta));
224 output->truth.phi.emplace_back(
static_cast<float>(
phi));
225 output->truth.charge.emplace_back(
static_cast<float>(
charge));
226 output->truth.extrapolatedStationMask.emplace_back(stationMask);
227 output->truth.m1Eta.emplace_back(positions[0].
eta);
228 output->truth.m1Phi.emplace_back(positions[0].
phi);
229 output->truth.m2Eta.emplace_back(positions[1].
eta);
230 output->truth.m2Phi.emplace_back(positions[1].
phi);
231 output->truth.m3Eta.emplace_back(positions[2].
eta);
232 output->truth.m3Phi.emplace_back(positions[2].
phi);
233 output->truth.matched.emplace_back(0U);
234 output->truth.matchedCandidateIndex.emplace_back(-1);
236 output->truth.unmatchedReason.emplace_back(
237 stationMask == 0x7U
238 ? static_cast<std::uint8_t>(
240 : static_cast<std::uint8_t>(
242 output->truth.wireSegmentMatched.emplace_back(0U);
243 output->truth.stripSegmentMatched.emplace_back(0U);
244 output->truth.matchedWireSegmentIndex.emplace_back(-1);
245 output->truth.matchedStripSegmentIndex.emplace_back(-1);
246 output->truth.wireSegmentMatchResidual.emplace_back(
248 output->truth.stripSegmentMatchResidual.emplace_back(
250 return StatusCode::SUCCESS;
251 };
252
254 if (event == nullptr) continue;
255#if __has_include("HepMC3/GenEvent.h")
256 for (
const auto& particle :
event->particles()) {
258 }
259#else
260 for (
auto particle =
event->particles_begin();
261 particle !=
event->particles_end(); ++particle) {
263 }
264#endif
265 }
266
267 for (const TgcL0Candidate& candidate : *candidates) {
268 output->candidates.subdetectorId.emplace_back(candidate.subdetectorId);
269 output->candidates.triggerSector.emplace_back(candidate.sectorId);
270 output->candidates.readoutSector.emplace_back(candidate.readoutSector);
271 output->candidates.bcTag.emplace_back(candidate.bcTag);
272 output->candidates.stationMask.emplace_back(candidate.positionStationMask);
273 output->candidates.wireStationMask.emplace_back(candidate.wireStationMask);
274 output->candidates.stripStationMask.emplace_back(candidate.stripStationMask);
275 output->candidates.eta.emplace_back(candidate.eta);
276 output->candidates.phi.emplace_back(candidate.phi);
277 output->candidates.deltaTheta.emplace_back(candidate.deltaTheta);
278 output->candidates.deltaPhi.emplace_back(candidate.deltaPhi);
279 output->candidates.truthIndex.emplace_back(-1);
280 }
281
282 std::vector<CandidateMatch> candidateMatches;
283 std::vector<bool> hasCandidateInWindow(
output->truth.pt.size(),
false);
284 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
285 if (
output->truth.extrapolatedStationMask[truth] != 0x7U)
continue;
286 const StationPositions positions = truthPositions(*output, truth);
287 for (std::size_t candidate = 0U; candidate <
candidates->size();
288 ++candidate) {
289 const TgcL0Candidate& inputCandidate = (*candidates)[candidate];
291 inputCandidate.bcTag !=
static_cast<std::uint16_t
>(
m_requiredBcTag.value())) {
292 continue;
293 }
294 if (
output->truth.eta[truth] * inputCandidate.eta < 0.F)
continue;
295 const float residual = meanDeltaR(inputCandidate, positions);
296 if (!std::isfinite(residual) || residual >
m_maxMeanDeltaR.value())
continue;
297 hasCandidateInWindow[truth] = true;
298 candidateMatches.push_back({
residual, truth, candidate});
299 }
300 }
302 [](const CandidateMatch& lhs, const CandidateMatch& rhs) {
303 if (lhs.residual != rhs.residual) {
304 return lhs.residual < rhs.residual;
305 }
306 if (lhs.truth != rhs.truth)
return lhs.truth < rhs.truth;
307 return lhs.candidate < rhs.candidate;
308 });
309 for (
const CandidateMatch&
match : candidateMatches) {
310 if (
output->truth.matchedCandidateIndex[
match.truth] >= 0 ||
311 output->candidates.truthIndex[
match.candidate] >= 0) {
312 continue;
313 }
316 static_cast<int>(
match.candidate);
321 static_cast<int>(
match.truth);
322 }
323 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
324 if (
output->truth.matched[truth] != 0U ||
325 output->truth.extrapolatedStationMask[truth] != 0x7U) {
326 continue;
327 }
328 output->truth.unmatchedReason[truth] =
static_cast<std::uint8_t
>(
329 hasCandidateInWindow[truth]
332 }
333
334 for (const TgcL0Segment& segment : *segments) {
335 output->segments.subdetectorId.emplace_back(segment.subdetectorId);
336 output->segments.triggerSector.emplace_back(segment.triggerSector);
337 output->segments.bcTag.emplace_back(segment.bcTag);
338 output->segments.projection.emplace_back(
339 static_cast<std::uint8_t>(segment.projection));
340 output->segments.stationMask.emplace_back(segment.stationMask);
341 output->segments.summedQuality.emplace_back(segment.summedQuality);
342 output->segments.nStations.emplace_back(segment.nStations);
343 output->segments.eta.emplace_back(segment.eta);
344 output->segments.phi.emplace_back(segment.phi);
345 output->segments.residual.emplace_back(segment.residual);
346 output->segments.outputResidual.emplace_back(segment.outputResidual);
347 output->segments.consistency.emplace_back(segment.consistency);
348 output->segments.pivotChannel.emplace_back(segment.pivotChannel);
349 output->segments.truthIndex.emplace_back(-1);
350 output->segments.truthMatchResidual.emplace_back(
352 }
353
355 const float maximumResidual) {
356 const auto projectionValue = static_cast<std::uint8_t>(projection);
357 std::vector<SegmentMatch> possibleMatches;
358 for (std::size_t truth = 0U; truth <
output->truth.pt.size(); ++truth) {
359 for (std::size_t segment = 0U;
360 segment <
output->segments.projection.size(); ++segment) {
361 if (
output->segments.projection[segment] != projectionValue)
continue;
363 output->segments.bcTag[segment] !=
365 continue;
366 }
367 if (
output->truth.eta[truth] *
output->segments.eta[segment] < 0.F) {
368 continue;
369 }
370 const int station = pivotStation(
output->segments.stationMask[segment]);
371 if (station < 0 ||
372 (
output->truth.extrapolatedStationMask[truth] &
373 stationBit(static_cast<std::size_t>(station))) == 0U) {
374 continue;
375 }
376 const std::array<float, 3>
truthEta{
378 output->truth.m3Eta[truth]};
379 const std::array<float, 3>
truthPhi{
381 output->truth.m3Phi[truth]};
384 ? std::abs(
output->segments.eta[segment] - truthEta[station])
387 if (std::isfinite(residual) && residual <= maximumResidual) {
388 possibleMatches.push_back({
residual, truth, segment});
389 }
390 }
391 }
393 [](const SegmentMatch& lhs, const SegmentMatch& rhs) {
394 if (lhs.residual != rhs.residual) {
395 return lhs.residual < rhs.residual;
396 }
397 if (lhs.truth != rhs.truth)
return lhs.truth < rhs.truth;
398 return lhs.segment < rhs.segment;
399 });
400 for (
const SegmentMatch&
match : possibleMatches) {
401 int& truthSegment =
403 ?
output->truth.matchedWireSegmentIndex[
match.truth]
404 :
output->truth.matchedStripSegmentIndex[
match.truth];
405 if (truthSegment >= 0 ||
407 continue;
408 }
409 truthSegment =
static_cast<int>(
match.segment);
411 static_cast<int>(
match.truth);
414 output->truth.wireSegmentMatched[
match.truth] = 1U;
416 } else {
417 output->truth.stripSegmentMatched[
match.truth] = 1U;
419 }
420 }
421 };
422
427
430 ATH_MSG_ERROR(
"Refusing to record inconsistent validation data: "
432 return StatusCode::FAILURE;
433 }
434
435 SG::WriteHandle<TgcL0ValidationEvent> outputHandle{
m_outputKey, ctx};
437 return StatusCode::SUCCESS;
438}
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.
double charge(const T &p)
Gaudi::Property< float > m_maxMeanDeltaR
Gaudi::Property< int > m_requiredTruthStatus
SG::WriteHandleKey< TgcL0ValidationEvent > m_outputKey
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
Gaudi::Property< double > m_validationMaxR
SG::ReadHandleKey< TgcL0SegmentContainer > m_segmentKey
Gaudi::Property< float > m_maxStripSegmentDeltaPhi
Gaudi::Property< double > m_minAbsEta
Gaudi::Property< int > m_requiredBcTag
Gaudi::Property< double > m_maxAbsEta
Gaudi::Property< float > m_maxWireSegmentDeltaEta
Gaudi::Property< double > m_validationPlaneToleranceZ
ToolHandle< Trk::IExtrapolator > m_extrapolator
SG::ReadHandleKey< McEventCollection > m_truthEventKey
virtual bool isValid() override final
Can the handle be successfully dereferenced?
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
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.
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[