19 #include "Acts/Utilities/RangeXD.hpp"
25 double yCross = (y0 + loc.
location().z() * tanTheta);
26 return (loc.
minY() < yCross && yCross < loc. maxY());
43 const std::array<double, 2>& chambEdges) {
45 if (chambEdges[0] <= seedEdges[0] && chambEdges[1] >= seedEdges[1]) {
49 else if (chambEdges[0]<= seedEdges[0]) {
50 return (chambEdges[1] - seedEdges[0]) / (seedEdges[1] - seedEdges[0]);
53 else if (chambEdges[1] >= seedEdges[1]) {
54 return (seedEdges[1] - chambEdges[0]) / (seedEdges[1] - seedEdges[0]);
57 else if (seedEdges[0] <= chambEdges[0] && seedEdges[1] >= chambEdges[1]) {
58 return (chambEdges[1] - chambEdges[0]) / (seedEdges[1] - seedEdges[0]);
69 return StatusCode::SUCCESS;
79 ATH_CHECK(writeMaxima.
record(std::make_unique<EtaHoughMaxContainer>()));
92 for (
auto& [station, stationHoughBuckets] :
data.houghSetups) {
94 for (
auto& bucket : stationHoughBuckets) {
98 writeMaxima->push_back(std::make_unique<HoughMaximum>(std::move(
max)));
102 std::stable_sort(writeMaxima->begin(), writeMaxima->end(),
104 return (*a->parentBucket()) < (*b->parentBucket());
106 return StatusCode::SUCCESS;
118 std::vector<HoughSetupForBucket>& buckets =
data.houghSetups[bucket->front()->msSector()];
121 Amg::Vector3D leftSide = globToLoc.translation() - (
hs.bucket->coveredMin() * Amg::Vector3D::UnitY());
122 Amg::Vector3D rightSide = globToLoc.translation() - (
hs.bucket->coveredMax() * Amg::Vector3D::UnitY());
126 for (
const std::shared_ptr<MuonR4::SpacePoint> & sp : *bucket) {
135 hs.searchWindowTanAngle = {tanThetaLeft, tanThetaRight};
140 for (
const std::shared_ptr<MuonR4::SpacePoint> & hit : *bucket){
142 double y0l = hit->localPosition().y() - hit->localPosition().z() * tanThetaLeft;
143 double y0r = hit->localPosition().y() - hit->localPosition().z() * tanThetaRight;
152 switch (hit->
type()){
162 return meas->
channelType() == sTgcIdHelper::sTgcChannelTypes::Strip;
170 HoughPlaneConfig
cfg;
177 data.houghPlane = std::make_unique<HoughPlane>(
cfg);
178 data.peakFinder = std::make_unique<ActsPeakFinderForMuon>(peakFinderCfg);
184 int expectedPrecisionChambers{0}, seenPrecisionChambers{0};
185 bool hasTrig =
false;
187 std::unordered_set<const MuonGMR4::MuonReadoutElement*> seenChambers{};
188 std::set<std::pair<int,int>> seenLayers;
189 using enum Acts::TrapezoidVolumeBounds::BoundValues;
192 std::array<double, 2> tubeExtend{halfX, -halfX};
194 auto addSeenHit = [&tubeExtend, &seenLayers, &seenChambers](
const SpacePoint& sp,
196 const double sensorL,
197 const int mL,
const int layer){
199 seenChambers.insert(
re);
200 tubeExtend[0] =
std::min(tubeExtend[0], sp.localPosition().x() - sensorL);
201 tubeExtend[1] =
std::max(tubeExtend[1], sp.localPosition().x() + sensorL);
203 for (
const SpacePoint* SP : maximum.hitIdentifiers){
204 ATH_MSG_VERBOSE(__func__<<
"() - "<<__LINE__<<
" Maximum has associated hit in "
211 addSeenHit(*SP,
re, 0.5*
re->activeTubeLength(dc->measurementHash()),
212 re->multilayer(), dc->tubeLayer());
216 addSeenHit(*SP,
re, 0.5*
re->stripLayer(clust->measurementHash()).design().stripLength(clust->channelNumber()),
217 re->multilayer(), clust->gasGap());
221 addSeenHit(*SP,
re, 0.5*
re->stripLayer(clust->measurementHash()).design().stripLength(clust->channelNumber()),
222 re->multilayer(), clust->gasGap());
232 ATH_MSG_VERBOSE(__func__<<
"() - "<<__LINE__<<
" maximum does not cross "
233 <<
m_idHelperSvc->toStringDetEl(muonChamber.readoutEle()->identify()));
237 <<
m_idHelperSvc->toStringDetEl(muonChamber.readoutEle()->identify())<<
", "
238 <<(*muonChamber.bounds())<<
", @: "<<
Amg::toString(muonChamber.location()));
243 const bool hasHit = seenChambers.count(muonChamber.readoutEle());
249 ++seenPrecisionChambers;
254 }
else if (precTech) {
257 const double lowL = muonChamber.width(maximum.y + muonChamber.location().z() * maximum.x);
258 const std::array<double, 2> chambEdges{muonChamber.location().x() - lowL,
259 muonChamber.location().x() + lowL};
262 if (coverage < 0.95){
263 ATH_MSG_VERBOSE(__func__<<
"() - "<<__LINE__<<
" Reject chamber due to partial coverage: "<<coverage
264 <<
", chamb extend: "<<chambEdges[0]<<
"-"<<chambEdges[1]
265 <<
", tube extend: "<<tubeExtend[0]<<
"-"<<tubeExtend[1]
266 <<
", "<<(*muonChamber.readoutEle()->msSector()->bounds()));
272 ++expectedPrecisionChambers;
276 int minLayers = seenLayers.size() / 2 + 1;
278 int minSeenPrecisionChambers = (expectedPrecisionChambers > 1) + 1;
282 minSeenPrecisionChambers = 1;
285 ", seen prec: "<<seenPrecisionChambers<<
", required: "<<minSeenPrecisionChambers
286 <<
" -- layers: "<<seenLayers.size()<<
", required: "<<
minLayers);
287 return seenPrecisionChambers >= minSeenPrecisionChambers && (
int)seenLayers.size() >=
minLayers;
318 data.currAxisRanges = Acts::HoughTransformUtils::HoughAxisRanges{
319 searchStartTanTheta, searchEndTanTheta, searchStart, searchEnd};
321 data.houghPlane->reset();
322 for (
const SpacePointBucket::value_type& hit : *(bucket.
bucket)) {
325 auto maxima =
data.peakFinder->findPeaks(*(
data.houghPlane),
data.currAxisRanges);
328 "#eta Hough accumulator");
330 if (maxima.empty()) {
332 <<
":\n Mean tanTheta was "<<tanThetaMean
333 <<
" and my intercept "<<chamberCenter
336 <<
". The bucket found a search range of ("
341 <<
") , and my final search range is ["
342 <<searchStartTanTheta<<
" - "<<searchEndTanTheta
343 <<
"] and ["<<searchStart<<
" - "<<searchEnd
351 std::set<HoughHitType> seenHits;
354 for (
const auto&
max : maxima) {
357 unsigned int nPrec{0};
358 auto toBins = [&
data](
double x,
double y){
359 return std::make_pair(
364 auto accumulatorBins = toBins(
max.x,
max.y);
365 for (
const HoughHitType& hit :
data.houghPlane->hitIds(accumulatorBins.first, accumulatorBins.second)) {
366 auto res = seenHits.emplace(hit);
377 std::vector<HoughHitType> hitList{
max.hitIdentifiers.begin(),
max.hitIdentifiers.end()};
394 switch (
chamber.readoutEle()->detectorType()) {
403 primitives.push_back(
MuonValR4::drawLabel(
std::format(
"Missed seed - score {}, layer score {}, comprising {} measurements ",
data.houghPlane->nHits(accumulatorBins.first, accumulatorBins.second),
data.houghPlane->nLayers(accumulatorBins.first, accumulatorBins.second),hitList.size()),0.05,0.03,12));
405 primitivesForAcc.push_back(
MuonValR4::drawLabel(
std::format(
"Missed seed - score {}, layer score {}, comprising {} measurements ",
data.houghPlane->nHits(accumulatorBins.first, accumulatorBins.second),
data.houghPlane->nLayers(accumulatorBins.first, accumulatorBins.second),hitList.size()),0.05,0.03,12));
407 "MissedAccumulator", std::move(primitivesForAcc));
408 m_visionTool->visualizeSeed(ctx, seed,
"Missed seed",std::move(primitives));
414 size_t nHits = hitList.size();
419 std::ranges::stable_sort(hitList,
sorter);
431 primitives.push_back(
MuonValR4::drawLabel(
std::format(
"score {}, layer score {}, comprising {} measurements. wx = {:.2f}, wy = {:.1f} ",
data.houghPlane->nHits(accumulatorBins.first, accumulatorBins.second),
data.houghPlane->nLayers(accumulatorBins.first, accumulatorBins.second),hitList.size(),
max.wx,
max.wy),0.05,0.03,12));
432 primitivesForAcc.push_back(
MuonValR4::drawLabel(
std::format(
"score {}, layer score {}, comprising {} measurements. wx = {:.2f}, wy = {:.1f} ",
data.houghPlane->nHits(accumulatorBins.first, accumulatorBins.second),
data.houghPlane->nLayers(accumulatorBins.first, accumulatorBins.second),hitList.size(),
max.wx /
m_targetResoTanTheta,
max.wy /
m_targetResoIntercept),0.05,0.03,12));
434 m_visionTool->visualizeAccumulator(ctx, *
data.houghPlane,
data.currAxisRanges, {max},
"#eta Hough accumulator", std::move(primitivesForAcc));
435 m_visionTool->visualizeSeed(ctx, seed,
"#eta-HoughSeed", std::move(primitives));
441 using namespace std::placeholders;
456 const unsigned precisionLayerIndex = (dc->readoutElement()->multilayer() * 10 + dc->tubeLayer());
471 for (
const SpacePointBucket::value_type& hit : *bucket.
bucket) {
472 if (!hit->measuresEta()) {
473 hitList.push_back(hit.get());