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 
14 
15 namespace {
16  constexpr double resetVal = 1.e10;
17 }
18 
19 
20 namespace MuonR4{
21 using namespace SegmentFit;
22 
24  ATH_CHECK(m_geoCtxKey.initialize());
25  ATH_CHECK(m_maxima.initialize());
26  ATH_CHECK(m_segmentSeeds.initialize());
27  ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
28  return StatusCode::SUCCESS;
29 }
31  HoughPlaneConfig cfg;
32  cfg.nBinsX = m_nBinsTanPhi;
33  cfg.nBinsY = m_nBinsIntercept;
34  // configure the peak finder for the phi-extension.
35  // Expect "shallow" maxima with 2-3 hits.
36  ActsPeakFinderForMuonCfg peakFinderCfg;
37  peakFinderCfg.fractionCutoff = 0.4;
38  peakFinderCfg.threshold = 2; // 2D spacepoints receive a weight of 2
39  peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
40  data.houghPlane = std::make_unique<HoughPlane>(cfg);
41  data.peakFinder = std::make_unique<ActsPeakFinderForMuon>(peakFinderCfg);
42 }
43 
44 int PhiHoughTransformAlg::countIncompatibleEtaHits(const ActsPeakFinderForMuon::Maximum& phiMaximum,
45  const HoughMaximum& etaMaximum) const {
46  std::unordered_map<const xAOD::UncalibratedMeasurement*, bool> foundEtas;
47  // loop over the original eta maximum and check all hits measuring the eta-coordinate
48  for (auto& hit : etaMaximum.getHitsInMax()) {
49  if (hit->measuresEta() && hit->measuresPhi()) {
50  auto [iter, added] = foundEtas.emplace(hit->primaryMeasurement(), false);
51  // test if the PRD for the eta-measurement appears on at least
52  // one space point of this phi extension.
53  // This will be done for all space points containing the eta-PRD
54  iter->second |= phiMaximum.hitIdentifiers.count(hit);
55  }
56  }
57  // count the number of eta PRD not compatible with this extension.
58  return std::ranges::count_if(foundEtas,
59  [](const std::pair<const xAOD::UncalibratedMeasurement*, bool>& p) {
60  return !p.second;
61  });
62 }
63 std::unique_ptr<SegmentSeed>
65  const ActsPeakFinderForMuon::Maximum & phiMax) const {
66  // book a new hit list
67  std::vector<HoughHitType> hitsOnMax{};
68  // copy the pure eta hits onto the hit list
69  std::ranges::copy_if(etaMax.getHitsInMax(),
70  std::back_inserter(hitsOnMax), [](const HoughHitType &hit){
71  return (hit->measuresEta() && !hit->measuresPhi());
72  });
73  // and then add all hits (2D and pure phi) from the phi-extension to it
74  hitsOnMax.insert(hitsOnMax.end(), phiMax.hitIdentifiers.begin(), phiMax.hitIdentifiers.end());
75  // use this to construct the segment seed
77  std::ranges::stable_sort(hitsOnMax, sorter);
78  return std::make_unique<SegmentSeed>(etaMax.tanTheta(), etaMax.interceptY(), phiMax.x, phiMax.y, hitsOnMax.size(), std::move(hitsOnMax), etaMax.parentBucket());
79 }
80 
82  const HoughMaximum & maximum,
83  HoughEventData& eventData) const{
84  // reset the event data
85  eventData.phiHitsOnMax = 0;
86  eventData.searchSpaceTanAngle = std::make_pair(resetVal, -resetVal);
87  eventData.searchSpaceIntercept = std::make_pair(resetVal, -resetVal);
88  // loop over the measurements on the maximum
89  for (auto hit : maximum.getHitsInMax()) {
90  // reject the pure eta measurements - not relevant here
91  if (!hit->measuresPhi()) {
92  continue;
93  }
94  // find the direction of the IP viewed from the sector frame
95  const Amg::Vector3D extrapDir = (hit->localPosition() - hit->msSector()->globalToLocalTrans(gctx).translation()).unit();
96  ATH_MSG_VERBOSE("Direction "<<Amg::toString(extrapDir));
97  // express the x location of our phi hits on the chamber plane (z = 0) when projecting from the beam spot
98  std::optional<double> dummyIntercept = Amg::intersect<3>(hit->localPosition(), extrapDir, Amg::Vector3D::UnitZ(),0);
99  double x0 = (hit->localPosition() + dummyIntercept.value_or(0) * extrapDir).x();
100  // now we can obtain the most likely tan(phi) via the pointing vector from the origin to our hit
101  double tanPhi = houghTanPhi(extrapDir);
102  // update our search space with this info
103  eventData.updateSearchWindow(eventData.searchSpaceTanAngle, tanPhi);
104  eventData.updateSearchWindow(eventData.searchSpaceIntercept, x0);
105  // and increment the hit counter
106  ++eventData.phiHitsOnMax;
107  }
108  // now use the results from the individual hits to define an axis range adapted to the binning we desire
109  double chamberCenter = (eventData.searchSpaceIntercept.second + eventData.searchSpaceIntercept.first) * 0.5;
110  double searchStart = chamberCenter - 0.5 * eventData.houghPlane->nBinsY() * m_targetResoIntercept;
111  double searchEnd = chamberCenter + 0.5 * eventData.houghPlane->nBinsY() * m_targetResoIntercept;
112  // Protection for very wide buckets - if the search space does not cover all of the bucket, widen the bin size
113  // so that we cover everything
114  searchStart = std::min(searchStart, eventData.searchSpaceIntercept.first- m_minSigmasSearchIntercept * m_targetResoIntercept);
115  searchEnd = std::max(searchEnd, eventData.searchSpaceIntercept.second + m_minSigmasSearchIntercept * m_targetResoIntercept);
116  // also treat tan(phi)
117  double tanPhiMean = 0.5 * (eventData.searchSpaceTanAngle.first + eventData.searchSpaceTanAngle.second);
118  double searchStartTanPhi = tanPhiMean - 0.5 * eventData.houghPlane->nBinsX() * m_targetResoTanPhi;
119  double searchEndTanPhi = tanPhiMean + 0.5* eventData.houghPlane->nBinsX() * m_targetResoTanPhi;
120  searchStartTanPhi = std::min(searchStartTanPhi, eventData.searchSpaceTanAngle.first- m_minSigmasSearchTanPhi * m_targetResoTanPhi);
121  searchEndTanPhi = std::max(searchEndTanPhi, eventData.searchSpaceTanAngle.second + m_minSigmasSearchTanPhi * m_targetResoTanPhi);
122 
123  // and update the axis ranges for the search space according to our results
124  eventData.currAxisRanges =
125  Acts::HoughTransformUtils::HoughAxisRanges{searchStartTanPhi, searchEndTanPhi, searchStart, searchEnd};
126  ATH_MSG_VERBOSE("Accumulator search window: tanPhi: ["<<searchStartTanPhi<<";"<<searchEndTanPhi<<"], x0: ["
127  <<searchStart<<";"<<searchEnd<<"]");
128 }
129 
130 std::vector<ActsPeakFinderForMuon::Maximum>
132  HoughEventData & eventData, const HoughMaximum & maximum) const{
133  std::unordered_map<int, std::vector<ActsPeakFinderForMuon::Maximum>> rankedSeeds;
134  using namespace std::placeholders;
135  // reset the accumulator
136  eventData.houghPlane->reset();
137  // fill the accumulator with the phi measurements
138  for (auto hit : maximum.getHitsInMax()){
139  if (!hit->measuresPhi()) {
140  ATH_MSG_VERBOSE("Hit "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<" does not have a phi measurement");
141  continue;
142  }
143  ATH_MSG_VERBOSE("Fill hit "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<", "<<Amg::toString(hit->localPosition()));
144  eventData.houghPlane->fill<HoughHitType>(
145  hit, eventData.currAxisRanges,
147  std::bind(HoughHelpers::Phi::houghWidthStrip, _1, _2, m_targetResoIntercept), hit, 0,
148  // up-weigh 2D spacepoints w.r.t 1D phi hits to prevent
149  // discarding measurements known to be compatible in eta
150  (hit->measuresEta() ? 2.0 : 1.0) / (m_downWeightMultiplePrd? hit->nPhiInstanceCounts() : 1)
151  );
152  }
153  // run the peak finder
154  auto foundMaxPhi = eventData.peakFinder->findPeaks(*eventData.houghPlane, eventData.currAxisRanges);
155  if (m_visionTool.isEnabled()) {
156  m_visionTool->visualizeAccumulator(ctx,*eventData.houghPlane, eventData.currAxisRanges, foundMaxPhi,
157  "#phi accumulator");
158  }
159  // now rank the found peaks by the number of eta-compatible PRDs they would discard.
160  for (const auto& solution : foundMaxPhi) {
161  // for each solution, count how many eta PRDs would not be compatible with this phi extension.
162  // rank by the lowest number of such "holes".
163  rankedSeeds[countIncompatibleEtaHits(solution, maximum)].push_back(solution);
164  }
165  // return only the best solution(s)
166  auto best = rankedSeeds.begin();
167  // and apply the maximum hole cut.
168  if (best != rankedSeeds.end() && best->first <= m_maxEtaHolesOnMax){
169  return best->second;
170  }
171  return {};
172 }
173 std::unique_ptr<SegmentSeed>
175  // recovers cases of a single phi hit assuming a straight
176  // line extrapolation from the beam line to the phi measurement
177  std::vector<HoughHitType> hits{maximum.getHitsInMax()};
179  std::ranges::stable_sort(hits, sorter);
180  return std::make_unique<SegmentSeed>(maximum.tanTheta(), maximum.interceptY(),
181  data.searchSpaceTanAngle.first,
182  data.searchSpaceIntercept.first,
183  maximum.getCounts(),
184  std::move(hits), maximum.parentBucket());
185 
186 }
187 
188 
189 StatusCode PhiHoughTransformAlg::execute(const EventContext& ctx) const {
190 
191  // read the inputs
192  const EtaHoughMaxContainer* maxima{nullptr};
193  ATH_CHECK(SG::get(maxima, m_maxima, ctx));
194 
195  const ActsGeometryContext* gctx{nullptr};
196  ATH_CHECK(SG::get(gctx, m_geoCtxKey, ctx));
197 
198  // book the event data object
199  HoughEventData eventData{};
200 
201  // prepare the accumulator
202  prepareHoughPlane(eventData);
203 
204  // prepare our output collection
205  SG::WriteHandle writeMaxima{m_segmentSeeds, ctx};
206  ATH_CHECK(writeMaxima.record(std::make_unique<SegmentSeedContainer>()));
207 
208  // loop over the previously found eta-maxima for each station
209  for (const HoughMaximum* max : *maxima) {
210  // for each maximum, pre-process
211  ATH_MSG_VERBOSE("Search extra phi hits on maximum "<<max->msSector()->identString()<<", tanTheta: "<<max->tanTheta()
212  <<", y0: "<<max->interceptY());
213  if (m_visionTool.isEnabled() && msgLvl(MSG::VERBOSE)) {
214  for (const auto& truth : m_visionTool->getLabeledSegments(max->getHitsInMax())) {
215  const Parameters truthPars = localSegmentPars(*truth);
216  ATH_MSG_VERBOSE("Truth parameters "<<toString(truthPars)<<", tanPhi: "
217  <<houghTanPhi(Amg::dirFromAngles(truthPars[Acts::toUnderlying(ParamDefs::phi)],
218  truthPars[Acts::toUnderlying(ParamDefs::theta)])));
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:64
MuonR4::PhiHoughTransformAlg::prepareHoughPlane
void prepareHoughPlane(HoughEventData &data) const
prepare the hough plane once per event.
Definition: PhiHoughTransformAlg.cxx:30
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:23
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
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
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
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
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:35
postInclude.sorter
sorter
Definition: postInclude.SortInput.py:23
MuonR4::ActsPeakFinderForMuonCfg
Acts::HoughTransformUtils::PeakFinders::IslandsAroundMaxConfig ActsPeakFinderForMuonCfg
Definition: MuonHoughDefs.h:31
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::houghTanPhi
double houghTanPhi(const Amg::Vector3D &v)
: Returns the hough tanPhi [x] / [z]
Definition: SegmentFitterEventData.cxx:19
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:44
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:38
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
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:60
MuonR4::HoughEventData_impl::currAxisRanges
Acts::HoughTransformUtils::HoughAxisRanges currAxisRanges
Definition: HoughEventData.h:61
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:131
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:24
PhiHoughTransformAlg.h
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:81
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:189
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:33
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:174
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::SegmentFit::Parameters
AmgVector(Acts::toUnderlying(ParamDefs::nPars)) Parameters
Definition: MuonHoughDefs.h:45
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::HoughEventData_impl
Templated event data class for the phase-2 muon hough transform.
Definition: HoughEventData.h:22