ATLAS Offline Software
PhiHoughTransformAlg.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "PhiHoughTransformAlg.h"
6 
8 
9 
15 
16 namespace {
17  constexpr double resetVal = 1.e10;
18 }
19 
20 
21 namespace MuonR4{
22 using namespace SegmentFit;
23 
25  ATH_CHECK(m_geoCtxKey.initialize());
26  ATH_CHECK(m_maxima.initialize());
27  ATH_CHECK(m_segmentSeeds.initialize());
28  ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
29  return StatusCode::SUCCESS;
30 }
32  HoughPlaneConfig cfg;
33  cfg.nBinsX = m_nBinsTanPhi;
34  cfg.nBinsY = m_nBinsIntercept;
35  // configure the peak finder for the phi-extension.
36  // Expect "shallow" maxima with 2-3 hits.
37  ActsPeakFinderForMuonCfg peakFinderCfg;
38  peakFinderCfg.fractionCutoff = 0.4;
39  peakFinderCfg.threshold = 2; // 2D spacepoints receive a weight of 2
40  peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
41  data.houghPlane = std::make_unique<HoughPlane>(cfg);
42  data.peakFinder = std::make_unique<ActsPeakFinderForMuon>(peakFinderCfg);
43 }
44 
45 int PhiHoughTransformAlg::countIncompatibleEtaHits(const ActsPeakFinderForMuon::Maximum& phiMaximum,
46  const HoughMaximum& etaMaximum) const {
47  std::unordered_map<const xAOD::UncalibratedMeasurement*, bool> foundEtas;
48  // loop over the original eta maximum and check all hits measuring the eta-coordinate
49  for (auto& hit : etaMaximum.getHitsInMax()) {
50  if (hit->measuresEta() && hit->measuresPhi()) {
51  auto [iter, added] = foundEtas.emplace(hit->primaryMeasurement(), false);
52  // test if the PRD for the eta-measurement appears on at least
53  // one space point of this phi extension.
54  // This will be done for all space points containing the eta-PRD
55  iter->second |= phiMaximum.hitIdentifiers.count(hit);
56  }
57  }
58  // count the number of eta PRD not compatible with this extension.
59  return std::ranges::count_if(foundEtas,
60  [](const std::pair<const xAOD::UncalibratedMeasurement*, bool>& p) {
61  return !p.second;
62  });
63 }
64 std::unique_ptr<SegmentSeed>
66  const ActsPeakFinderForMuon::Maximum & phiMax) const {
67  // book a new hit list
68  std::vector<HoughHitType> hitsOnMax{};
69  // copy the pure eta hits onto the hit list
70  std::ranges::copy_if(etaMax.getHitsInMax(),
71  std::back_inserter(hitsOnMax), [](const HoughHitType &hit){
72  return (hit->measuresEta() && !hit->measuresPhi());
73  });
74  // and then add all hits (2D and pure phi) from the phi-extension to it
75  hitsOnMax.insert(hitsOnMax.end(), phiMax.hitIdentifiers.begin(), phiMax.hitIdentifiers.end());
76  // use this to construct the segment seed
78  std::ranges::stable_sort(hitsOnMax, sorter);
79  return std::make_unique<SegmentSeed>(etaMax.tanTheta(), etaMax.interceptY(), phiMax.x, phiMax.y, hitsOnMax.size(), std::move(hitsOnMax), etaMax.parentBucket());
80 }
81 
83  const HoughMaximum & maximum,
84  HoughEventData& eventData) const{
85  // reset the event data
86  eventData.phiHitsOnMax = 0;
87  eventData.searchSpaceTanAngle = std::make_pair(resetVal, -resetVal);
88  eventData.searchSpaceIntercept = std::make_pair(resetVal, -resetVal);
89  // loop over the measurements on the maximum
90  for (auto hit : maximum.getHitsInMax()) {
91  // reject the pure eta measurements - not relevant here
92  if (!hit->measuresPhi()) {
93  continue;
94  }
95  // find the direction of the IP viewed from the sector frame
96  const Amg::Vector3D extrapDir = (hit->positionInChamber() - hit->msSector()->globalToLocalTrans(gctx).translation()).unit();
97  ATH_MSG_VERBOSE("Direction "<<Amg::toString(extrapDir));
98  // express the x location of our phi hits on the chamber plane (z = 0) when projecting from the beam spot
99  std::optional<double> dummyIntercept = Amg::intersect<3>(hit->positionInChamber(), extrapDir, Amg::Vector3D::UnitZ(),0);
100  double x0 = (hit->positionInChamber() + dummyIntercept.value_or(0) * extrapDir).x();
101  // now we can obtain the most likely tan(phi) via the pointing vector from the origin to our hit
102  double tanPhi = houghTanPhi(extrapDir);
103  // update our search space with this info
104  eventData.updateSearchWindow(eventData.searchSpaceTanAngle, tanPhi);
105  eventData.updateSearchWindow(eventData.searchSpaceIntercept, x0);
106  // and increment the hit counter
107  ++eventData.phiHitsOnMax;
108  }
109  // now use the results from the individual hits to define an axis range adapted to the binning we desire
110  double chamberCenter = (eventData.searchSpaceIntercept.second + eventData.searchSpaceIntercept.first) * 0.5;
111  double searchStart = chamberCenter - 0.5 * eventData.houghPlane->nBinsY() * m_targetResoIntercept;
112  double searchEnd = chamberCenter + 0.5 * eventData.houghPlane->nBinsY() * m_targetResoIntercept;
113  // Protection for very wide buckets - if the search space does not cover all of the bucket, widen the bin size
114  // so that we cover everything
115  searchStart = std::min(searchStart, eventData.searchSpaceIntercept.first- m_minSigmasSearchIntercept * m_targetResoIntercept);
116  searchEnd = std::max(searchEnd, eventData.searchSpaceIntercept.second + m_minSigmasSearchIntercept * m_targetResoIntercept);
117  // also treat tan(phi)
118  double tanPhiMean = 0.5 * (eventData.searchSpaceTanAngle.first + eventData.searchSpaceTanAngle.second);
119  double searchStartTanPhi = tanPhiMean - 0.5 * eventData.houghPlane->nBinsX() * m_targetResoTanPhi;
120  double searchEndTanPhi = tanPhiMean + 0.5* eventData.houghPlane->nBinsX() * m_targetResoTanPhi;
121  searchStartTanPhi = std::min(searchStartTanPhi, eventData.searchSpaceTanAngle.first- m_minSigmasSearchTanPhi * m_targetResoTanPhi);
122  searchEndTanPhi = std::max(searchEndTanPhi, eventData.searchSpaceTanAngle.second + m_minSigmasSearchTanPhi * m_targetResoTanPhi);
123 
124  // and update the axis ranges for the search space according to our results
125  eventData.currAxisRanges =
126  Acts::HoughTransformUtils::HoughAxisRanges{searchStartTanPhi, searchEndTanPhi, searchStart, searchEnd};
127  ATH_MSG_VERBOSE("Accumulator search window: tanPhi: ["<<searchStartTanPhi<<";"<<searchEndTanPhi<<"], x0: ["
128  <<searchStart<<";"<<searchEnd<<"]");
129 }
130 
131 std::vector<ActsPeakFinderForMuon::Maximum>
133  HoughEventData & eventData, const HoughMaximum & maximum) const{
134  std::unordered_map<int, std::vector<ActsPeakFinderForMuon::Maximum>> rankedSeeds;
135  using namespace std::placeholders;
136  // reset the accumulator
137  eventData.houghPlane->reset();
138  // fill the accumulator with the phi measurements
139  for (auto hit : maximum.getHitsInMax()){
140  if (!hit->measuresPhi()) {
141  ATH_MSG_VERBOSE("Hit "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<" does not have a phi measurement");
142  continue;
143  }
144  ATH_MSG_VERBOSE("Fill hit "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<", "<<Amg::toString(hit->positionInChamber()));
145  eventData.houghPlane->fill<HoughHitType>(
146  hit, eventData.currAxisRanges,
148  std::bind(HoughHelpers::Phi::houghWidthStrip, _1, _2, m_targetResoIntercept), hit, 0,
149  // up-weigh 2D spacepoints w.r.t 1D phi hits to prevent
150  // discarding measurements known to be compatible in eta
151  (hit->measuresEta() ? 2.0 : 1.0) / (m_downWeightMultiplePrd? hit->nPhiInstanceCounts() : 1)
152  );
153  }
154  // run the peak finder
155  auto foundMaxPhi = eventData.peakFinder->findPeaks(*eventData.houghPlane, eventData.currAxisRanges);
156  if (m_visionTool.isEnabled()) {
157  m_visionTool->visualizeAccumulator(ctx,*eventData.houghPlane, eventData.currAxisRanges, foundMaxPhi,
158  "#phi accumulator");
159  }
160  // now rank the found peaks by the number of eta-compatible PRDs they would discard.
161  for (const auto& solution : foundMaxPhi) {
162  // for each solution, count how many eta PRDs would not be compatible with this phi extension.
163  // rank by the lowest number of such "holes".
164  rankedSeeds[countIncompatibleEtaHits(solution, maximum)].push_back(solution);
165  }
166  // return only the best solution(s)
167  auto best = rankedSeeds.begin();
168  // and apply the maximum hole cut.
169  if (best != rankedSeeds.end() && best->first <= m_maxEtaHolesOnMax){
170  return best->second;
171  }
172  return {};
173 }
174 std::unique_ptr<SegmentSeed>
176  // recovers cases of a single phi hit assuming a straight
177  // line extrapolation from the beam line to the phi measurement
178  std::vector<HoughHitType> hits{maximum.getHitsInMax()};
180  std::ranges::stable_sort(hits, sorter);
181  return std::make_unique<SegmentSeed>(maximum.tanTheta(), maximum.interceptY(),
182  data.searchSpaceTanAngle.first,
183  data.searchSpaceIntercept.first,
184  maximum.getCounts(),
185  std::move(hits), maximum.parentBucket());
186 
187 }
188 
189 
190 StatusCode PhiHoughTransformAlg::execute(const EventContext& ctx) const {
191 
192  // read the inputs
193  const EtaHoughMaxContainer* maxima{nullptr};
194  ATH_CHECK(SG::get(maxima, m_maxima, ctx));
195 
196  const ActsGeometryContext* gctx{nullptr};
197  ATH_CHECK(SG::get(gctx, m_geoCtxKey, ctx));
198 
199  // book the event data object
200  HoughEventData eventData{};
201 
202  // prepare the accumulator
203  prepareHoughPlane(eventData);
204 
205  // prepare our output collection
206  SG::WriteHandle writeMaxima{m_segmentSeeds, ctx};
207  ATH_CHECK(writeMaxima.record(std::make_unique<SegmentSeedContainer>()));
208 
209  // loop over the previously found eta-maxima for each station
210  for (const HoughMaximum* max : *maxima) {
211  // for each maximum, pre-process
212  ATH_MSG_VERBOSE("Search extra phi hits on maximum "<<max->msSector()->identString()<<", tanTheta: "<<max->tanTheta()
213  <<", y0: "<<max->interceptY());
214  if (m_visionTool.isEnabled() && msgLvl(MSG::VERBOSE)) {
215  for (const auto& truth : m_visionTool->getLabeledSegments(max->getHitsInMax())) {
216  const Parameters truthPars = localSegmentPars(*truth);
217  ATH_MSG_VERBOSE("Truth parameters "<<toString(truthPars)<<", tanPhi: "
219  }
220  }
221  preProcessMaximum(*gctx, *max, eventData);
222  bool foundSolution=false;
223  // if we have enough hits, run a phi transform
224  if (eventData.phiHitsOnMax > 1){
225  std::vector<ActsPeakFinderForMuon::Maximum> rankedSeeds = findRankedSegmentSeeds(ctx, eventData, *max);
226  for (auto & phiSolution : rankedSeeds){
227  foundSolution = true;
228  const SegmentSeed* seed {writeMaxima->push_back(buildSegmentSeed(*max, phiSolution))};
229  if (m_visionTool.isEnabled()) {
230  m_visionTool->visualizeSeed(ctx, *seed, "#phi pattern seed");
231  }
232  }
233  }
234  ATH_MSG_VERBOSE("Solution found: "<<foundSolution);
235  // if we do not have at least two phi-hits for a proper transform:
236  if (!foundSolution){
237  // if we have a single phi hit, we can approximate the phi
238  // solution using the beam spot (as the IP is far).
239  // This is steered by a flag, and not appropriate for splashes
240  // or cosmics.
241  if (m_recoverSinglePhiWithBS && eventData.phiHitsOnMax == 1){
242  const SegmentSeed* singleMax{writeMaxima->push_back(recoverSinglePhiMax(eventData,*max))};
243  if (m_visionTool.isEnabled()) {
244  m_visionTool->visualizeSeed(ctx, *singleMax, "Single #phi hit recovery");
245  }
246  }
247  // otherwise we have no phi-solution, and we fall back to writing a 1D eta-maximum
248  else{
249  writeMaxima->push_back(std::make_unique<SegmentSeed>(*max));
250  }
251  }
252  }
253  // add the maxima for this station to the output
254 
255  return StatusCode::SUCCESS;
256 }
257 }
MuonR4::PhiHoughTransformAlg::buildSegmentSeed
std::unique_ptr< SegmentSeed > buildSegmentSeed(const HoughMaximum &etaMax, const MuonR4::ActsPeakFinderForMuon::Maximum &phiMax) const
constructs a segment seed from an eta maximum and a phi-extension.
Definition: PhiHoughTransformAlg.cxx:65
MuonR4::PhiHoughTransformAlg::prepareHoughPlane
void prepareHoughPlane(HoughEventData &data) const
prepare the hough plane once per event.
Definition: PhiHoughTransformAlg.cxx:31
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
createLinkingScheme.iter
iter
Definition: createLinkingScheme.py:62
MuonR4::PhiHoughTransformAlg::initialize
virtual StatusCode initialize() override
Definition: PhiHoughTransformAlg.cxx:24
MuonR4::HoughMaximum::tanTheta
double tanTheta() const
getter
Definition: HoughMaximum.h:36
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
MuonR4::HoughEventData_impl::peakFinder
std::unique_ptr< peakFinder_t > peakFinder
Definition: HoughEventData.h:37
MuonR4::SegmentFit::Parameters
AmgVector(toInt(ParamDefs::nPars)) Parameters
Definition: MuonHoughDefs.h:48
MuonR4::HoughEventData_impl::phiHitsOnMax
size_t phiHitsOnMax
Definition: HoughEventData.h:67
MuonR4::HoughEventData_impl::searchSpaceTanAngle
std::pair< double, double > searchSpaceTanAngle
Definition: HoughEventData.h:65
MuonR4::SpacePointPerLayerSorter
The SpacePointPerLayerSorter sort two given space points by their layer Identifier.
Definition: SpacePointPerLayerSorter.h:15
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
EventPrimitivesHelpers.h
HoughHelperFunctions.h
xAOD::etaMax
etaMax
Definition: HIEventShape_v2.cxx:46
MuonR4::HoughHelpers::Phi::houghParamStrip
double houghParamStrip(double tanPhi, const MuonR4::HoughHitType &dc)
straight line parametrisation for strip detector measurements, in the x-direction
Definition: HoughHelperFunctions.cxx:31
postInclude.sorter
sorter
Definition: postInclude.SortInput.py:23
MuonR4::ActsPeakFinderForMuonCfg
Acts::HoughTransformUtils::PeakFinders::IslandsAroundMaxConfig ActsPeakFinderForMuonCfg
Definition: MuonHoughDefs.h:26
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
SpectrometerSector.h
x
#define x
MuonR4::HoughMaximum::getHitsInMax
const std::vector< HitType > & getHitsInMax() const
getter
Definition: HoughMaximum.h:46
MuonR4::SegmentFit::ParamDefs::phi
@ phi
MuonR4::houghTanPhi
double houghTanPhi(const Amg::Vector3D &v)
: Returns the hough tanPhi [x] / [z]
Definition: SegmentFitterEventData.cxx:18
MuonR4::PhiHoughTransformAlg::countIncompatibleEtaHits
int countIncompatibleEtaHits(const MuonR4::ActsPeakFinderForMuon::Maximum &phiMaximum, const MuonR4::HoughMaximum &etaMaximum) const
helper to count the number of eta measurements that would be discarded for a given phi extension cand...
Definition: PhiHoughTransformAlg.cxx:45
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
MuonR4::HoughHelpers::Phi::houghWidthStrip
double houghWidthStrip(double tanPhi, const MuonR4::HoughHitType &dc, double targetReso)
Uncertainty parametrisation for strip measurements.
Definition: HoughHelperFunctions.cxx:34
MuonR4::HoughEventData_impl::searchSpaceIntercept
std::pair< double, double > searchSpaceIntercept
Definition: HoughEventData.h:63
SG::get
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.
Definition: ReadCondHandle.h:287
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
UtilFunctions.h
MuonR4::HoughMaximum::parentBucket
const SpacePointBucket * parentBucket() const
getter
Definition: HoughMaximum.h:51
python.utils.best
def best(iterable, priorities=[3, 2, 1, -1, 0])
Definition: DataQuality/DQUtils/python/utils.py:49
AnalysisUtils::copy_if
Out copy_if(In first, const In &last, Out res, const Pred &p)
Definition: IFilterUtils.h:30
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
MuonR4::SegmentFit::toString
std::string toString(const Parameters &pars)
Definition: SegmentFitterEventData.cxx:59
MuonR4::HoughEventData_impl::currAxisRanges
Acts::HoughTransformUtils::HoughAxisRanges currAxisRanges
Definition: HoughEventData.h:61
MuonR4::SegmentFit::ParamDefs::x0
@ x0
MuonR4::PhiHoughTransformAlg::findRankedSegmentSeeds
std::vector< MuonR4::ActsPeakFinderForMuon::Maximum > findRankedSegmentSeeds(const EventContext &ctx, HoughEventData &data, const MuonR4::HoughMaximum &maximum) const
perform a hough search for the most promising phi extension of an eta-maximum Performs a local hough ...
Definition: PhiHoughTransformAlg.cxx:132
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
ActsGeometryContext
Include the GeoPrimitives which need to be put first.
Definition: ActsGeometryContext.h:27
MuonR4::SpacePoint
The muon space point is the combination of two uncalibrated measurements one of them measures the eta...
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/MuonSpacePoint/SpacePoint.h:19
PhiHoughTransformAlg.h
MuonR4::SegmentFit::toInt
constexpr int toInt(const ParamDefs p)
Definition: MuonHoughDefs.h:42
MuonR4::HoughMaximum::getCounts
double getCounts() const
getter
Definition: HoughMaximum.h:43
WriteCaloSwCorrections.cfg
cfg
Definition: WriteCaloSwCorrections.py:23
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MuonR4::PhiHoughTransformAlg::preProcessMaximum
void preProcessMaximum(const ActsGeometryContext &gctx, const MuonR4::HoughMaximum &maximum, HoughEventData &data) const
pre-processing for a given input eta-maximum Counts potential phi-hits and defines the search space
Definition: PhiHoughTransformAlg.cxx:82
Amg::dirFromAngles
Amg::Vector3D dirFromAngles(const double phi, const double theta)
Constructs a direction vector from the azimuthal & polar angles.
Definition: GeoPrimitivesHelpers.h:299
MuonR4
This header ties the generic definitions in this package.
Definition: HoughEventData.h:16
MuonR4::PhiHoughTransformAlg::execute
virtual StatusCode execute(const EventContext &ctx) const override
Definition: PhiHoughTransformAlg.cxx:190
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:73
MuonR4::HoughMaximum
Data class to represent an eta maximum in hough space.
Definition: HoughMaximum.h:14
SegmentFitterEventData.h
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:21
MuonR4::SegmentSeed
Representation of a segment seed (a fully processed hough maximum) produced by the hough transform.
Definition: SegmentSeed.h:14
MuonR4::HoughMaximum::interceptY
double interceptY() const
getter
Definition: HoughMaximum.h:40
MuonR4::SegmentFit::localSegmentPars
Parameters localSegmentPars(const xAOD::MuonSegment &seg)
Returns the localSegPars decoration from a xAODMuon::Segment.
Definition: SegmentFitterEventData.cxx:32
MuonR4::PhiHoughTransformAlg::recoverSinglePhiMax
std::unique_ptr< SegmentSeed > recoverSinglePhiMax(HoughEventData &data, const MuonR4::HoughMaximum &maximum) const
extend an eta maximum with just a single attached phi measurement.
Definition: PhiHoughTransformAlg.cxx:175
MuonR4::HoughEventData_impl::houghPlane
std::unique_ptr< HoughPlane > houghPlane
Definition: HoughEventData.h:35
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:13
MuonR4::HoughEventData_impl::updateSearchWindow
static void updateSearchWindow(std::pair< double, double > &searchWindow, double value)
Updates a search space window to account for a value.
Definition: HoughEventData.h:29
SpacePointPerLayerSorter.h
MuonR4::SegmentFit::ParamDefs::theta
@ theta
MuonR4::HoughEventData_impl
Templated event data class for the phase-2 muon hough transform.
Definition: HoughEventData.h:22