ATLAS Offline Software
MuonLayerHoughTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include "AtlasHepMC/GenEvent.h"
8 #include "CxxUtils/sincos.h"
9 #include "GaudiKernel/ConcurrencyFlags.h"
20 namespace Muon {
21  using namespace MuonStationIndex;
22 
24  ATH_CHECK(m_idHelperSvc.retrieve());
25  m_ntechnologies = m_idHelperSvc->mdtIdHelper().technologyNameIndexMax() + 1;
26  ATH_CHECK(m_printer.retrieve());
27  ATH_CHECK(m_muonManagerKey.initialize());
28  ATH_CHECK(m_truthNames.initialize());
29  constexpr int nSelect = toInt(ChIndex::ChIndexMax);
30  // initialize cuts, if only one cut, use make_pair to avoid compiler issues, format is (position, cut)
31  m_selectors.resize(nSelect);
32  m_selectors[toInt(ChIndex::BIS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 5.9)}); // old values: 6.9; optimized: 7.9
33  m_selectors[toInt(ChIndex::BIL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 5.9)}); // old values: 6.9; optimized: 7.9
34  m_selectors[toInt(ChIndex::BMS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 7.9; optimized: 7.9
35  m_selectors[toInt(ChIndex::BML)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 7.9; optimized: 7.9
36  m_selectors[toInt(ChIndex::BOS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
37  m_selectors[toInt(ChIndex::BOL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
38  m_selectors[toInt(ChIndex::BEE)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 3.9)}); // old values: 5.9; optimized: 5.9
39  m_selectors[toInt(ChIndex::EIS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 5.9)}); // old values: 6.9; optimized: 7.9
40  m_selectors[toInt(ChIndex::EIL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 5.9)}); // old values: 6.9; optimized: 7.9
41  m_selectors[toInt(ChIndex::EMS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 7.9; optimized: 5.9
42  m_selectors[toInt(ChIndex::EML)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 7.9; optimized: 5.9
43  m_selectors[toInt(ChIndex::EOS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
44  m_selectors[toInt(ChIndex::EOL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
45  m_selectors[toInt(ChIndex::EES)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
46  m_selectors[toInt(ChIndex::EEL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 4.9)}); // old values: 4.9; optimized: 5.9
47 
48  m_selectorsLoose.resize(nSelect);
49  m_selectorsLoose[toInt(ChIndex::BIS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 2.9; optimized: 3.9
50  m_selectorsLoose[toInt(ChIndex::BIL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 3.9
51  m_selectorsLoose[toInt(ChIndex::BMS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 4.9; optimized: 2.9
52  m_selectorsLoose[toInt(ChIndex::BML)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 4.9; optimized: 2.9
53  m_selectorsLoose[toInt(ChIndex::BOS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 2.9; optimized: 2.9
54  m_selectorsLoose[toInt(ChIndex::BOL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 2.9
55  m_selectorsLoose[toInt(ChIndex::BEE)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 3.9; optimized: 2.9
56  m_selectorsLoose[toInt(ChIndex::EIS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 4.9; optimized: 3.9
57  m_selectorsLoose[toInt(ChIndex::EIL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 1.9)}); // old values: 4.9; optimized: 3.9
58  m_selectorsLoose[toInt(ChIndex::EMS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 5.9; optimized: 2.9
59  m_selectorsLoose[toInt(ChIndex::EML)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 5.9; optimized: 2.9
60  m_selectorsLoose[toInt(ChIndex::EOS)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 2.9
61  m_selectorsLoose[toInt(ChIndex::EOL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 2.9
62  m_selectorsLoose[toInt(ChIndex::EES)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 2.9
63  m_selectorsLoose[toInt(ChIndex::EEL)] = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 2.9)}); // old values: 2.9; optimized: 2.9
64 
65  return StatusCode::SUCCESS;
66  }
67 
68  std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> MuonLayerHoughTool::find(
69  const std::vector<const MdtPrepDataCollection*>& mdtCols, const std::vector<const CscPrepDataCollection*>& cscCols,
70  const std::vector<const TgcPrepDataCollection*>& tgcCols, const std::vector<const RpcPrepDataCollection*>& rpcCols,
71  const MuonSegmentCombinationCollection*, const EventContext& ctx) const {
72 
73  using namespace MuonStationIndex;
74  initializeSectorMapping(ctx);
75  State state;
76  ATH_MSG_DEBUG("MuonLayerHoughTool::find");
77 
78  // create structure to hold data per sector and set the sector indices
79  state.houghDataPerSectorVec->vec.resize(16);
80  for (unsigned int i = 0; i < state.houghDataPerSectorVec->vec.size(); ++i) state.houghDataPerSectorVec->vec[i].sector = i + 1;
81 
82  // return DetectorRegionIndex and sectorLayerHash
83  auto getHashes = [this](const Identifier& id) {
84  DetRegIdx regionIndex = m_idHelperSvc->regionIndex(id);
85  MuonStationIndex::LayerIndex layerIndex = m_idHelperSvc->layerIndex(id);
86  unsigned int sectorLayerHash = MuonStationIndex::sectorLayerHash(regionIndex, layerIndex);
87  return std::make_pair(regionIndex, sectorLayerHash);
88  };
89 
90  for (const MdtPrepDataCollection* col : mdtCols) {
91  if (!col) continue;
92  Identifier id = col->identify();
93  int sector = m_idHelperSvc->sector(id);
94  auto hashes = getHashes(id);
95  fill(ctx, state.truthHits, *col, state.houghDataPerSectorVec->vec[sector - 1].hitVec[hashes.second]);
96  }
97 
98  for (const RpcPrepDataCollection* col : rpcCols) {
99  if (!col) continue;
100  Identifier id = col->identify();
101  int sector = m_idHelperSvc->sector(id);
102  auto hashes = getHashes(id);
103  fill(ctx, state.truthHits, *col, state.houghDataPerSectorVec->vec[sector - 1].hitVec[hashes.second],
104  state.houghDataPerSectorVec->vec[sector - 1].phiHitVec[toInt(hashes.first)]);
105  }
106  for (const CscPrepDataCollection* col : cscCols) {
107  if (!col) continue;
108  const Identifier id = col->identify();
109  int sector = m_idHelperSvc->sector(id);
110  auto hashes = getHashes(id);
111  fill(ctx, state.truthHits, *col, state.houghDataPerSectorVec->vec[sector - 1].hitVec[hashes.second],
112  state.houghDataPerSectorVec->vec[sector - 1].phiHitVec[toInt(hashes.first)]);
113  }
114  auto hashInSector = [this](IdentifierHash hash, int sector, unsigned int sectorLayerHash) {
115  constexpr int tgcTech = toInt(TechnologyIndex::TGC);
116  const std::vector<IdentifierHash>& hashes = m_collectionsPerSector[sector - 1].technologyRegionHashVecs[tgcTech][sectorLayerHash];
117  return std::binary_search(hashes.begin(), hashes.end(), hash);
118  };
119 
120  for (const TgcPrepDataCollection* col : tgcCols) {
121  if (!col) continue;
122  Identifier id = col->identify();
123  int sector = m_idHelperSvc->sector(id);
124  auto hashes = getHashes(id);
125  // fill current sector
126  fill(ctx, state.truthHits, state.houghDataPerSectorVec->tgcClusteringObjs, *col,
127  state.houghDataPerSectorVec->vec[sector - 1].hitVec[hashes.second],
128  state.houghDataPerSectorVec->vec[sector - 1].phiHitVec[toInt(hashes.first)], sector);
129 
130  // fill neighbours if in overlap
131  int neighbourSectorDown = sector == 1 ? 16 : sector - 1;
132  if (hashInSector(col->identifyHash(), neighbourSectorDown, hashes.second))
133  fill(ctx, state.truthHits, state.houghDataPerSectorVec->tgcClusteringObjs, *col,
134  state.houghDataPerSectorVec->vec[neighbourSectorDown - 1].hitVec[hashes.second],
135  state.houghDataPerSectorVec->vec[neighbourSectorDown - 1].phiHitVec[toInt(hashes.first)], neighbourSectorDown);
136 
137  int neighbourSectorUp = sector == 16 ? 1 : sector + 1;
138  if (hashInSector(col->identifyHash(), neighbourSectorUp, hashes.second))
139  fill(ctx, state.truthHits, state.houghDataPerSectorVec->tgcClusteringObjs, *col,
140  state.houghDataPerSectorVec->vec[neighbourSectorUp - 1].hitVec[hashes.second],
141  state.houghDataPerSectorVec->vec[neighbourSectorUp - 1].phiHitVec[toInt(hashes.first)], neighbourSectorUp);
142  }
143 
144  return analyse(state);
145  }
146 
147  std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> MuonLayerHoughTool::find(
148  const MdtPrepDataContainer* mdtCont, const CscPrepDataContainer* cscCont, const TgcPrepDataContainer* tgcCont,
149  const RpcPrepDataContainer* rpcCont, const sTgcPrepDataContainer* stgcCont, const MMPrepDataContainer* mmCont,
150  const EventContext& ctx) const {
151  initializeSectorMapping(ctx);
152  State state;
153  ATH_MSG_DEBUG("MuonLayerHoughTool::analyse");
154 
155  state.houghDataPerSectorVec->vec.resize(16);
156 
157  // loops over all sectors, contains hashes for technology and region and chamber (?)
158  for (const CollectionsPerSector& sit : m_collectionsPerSector) {
159  ATH_MSG_DEBUG("analyse: Filling hits sector " << sit.sector);
160 
161  HoughDataPerSector& houghData = state.houghDataPerSectorVec->vec[sit.sector - 1];
162  houghData.sector = sit.sector;
163 
164  // fill hits for this sector -> hitsVec and PhiHitsVec are known now
165  fillHitsPerSector(ctx, state, sit.sector, sit, mdtCont, cscCont, tgcCont, rpcCont, stgcCont, mmCont);
166  }
167  return analyse(state);
168  }
169 
170  std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> MuonLayerHoughTool::analyse(
171  State& state) const {
172  auto patternCombis = std::make_unique<MuonPatternCombinationCollection>();
173 
174  // loop over data and fill the hough transform
175  for (auto& houghData : state.houghDataPerSectorVec->vec) {
176  ATH_MSG_DEBUG("analyse: Filling Hough sector " << houghData.sector);
177 
178 
179  // loop over all possible station layers in the sector and run the eta transform
180  for (unsigned int layerHash = 0; layerHash < MuonStationIndex::sectorLayerHashMax(); ++layerHash) {
181  using namespace MuonStationIndex;
182  // get hits for layer, skip empty layers
183  HitVec& hits = houghData.hitVec.at(layerHash);
184  if (hits.empty()) continue;
185 
186  // decompose hash, calculate indices etc
187  const auto [region , layer] = decomposeSectorLayerHash(layerHash);
188  StIndex index = toStationIndex(region, layer);
189 
190  // get Hough transform
192  state.houghDataPerSectorVec->detectorHoughTransforms.hough(houghData.sector, region, layer);
193 
194  ATH_MSG_VERBOSE("analyse: Filling Summary: loc s" << houghData.sector << " " << regionName(region) << " "
195  << layerName(layer)<<" -> hash: "<<layerHash << " -> stIndex: "
196  << stName(index) << " etaHits: " << hits.size());
197 
198  // look for maxima using hough in eta per layer
199  if (!findMaxima(state.seedMaxima, hough, hits, houghData.maxVec.at(layerHash)) ||
200  houghData.maxVec.at(layerHash).empty())
201  continue;
202 
203  ++houghData.nlayersWithMaxima[toInt(region)];
204  houghData.nmaxHitsInRegion[toInt(region)] += houghData.maxVec[layerHash].front()->max;
205 
206  ATH_MSG_VERBOSE("analyse: Eta maxima Summary: loc s"
207  << houghData.sector << " " << regionName(region) << " "
208  << layerName(layer) << " -> stIndex: " << stName(index)
209  << " hash: " << layerHash << " nMaxima: " << houghData.maxVec[layerHash].size());
210  } // loop over layerHash -> maxima per layer in eta are known now
211  } // loop over sectors
212 
213  if (m_useSeeds) {
214  std::vector<Road> roads;
215  buildRoads(state.seedMaxima, state.houghDataPerSectorVec->detectorHoughTransforms,
216  state.houghDataPerSectorVec, roads);
217 
218  // create association map
219  ATH_MSG_DEBUG("analyse: Building pattern combinations using roads " << roads.size());
220  for (auto& road : roads) {
221  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MuonLayerHoughTool::MaximumVec> phiEtaAssMap;
222  MuonLayerHoughTool::RegionMaximumVec unassociatedEtaMaxima;
223 
224  int sector = road.seed->hough->m_descriptor.sector;
225  using namespace MuonStationIndex;
226  ChIndex chIndex = road.seed->hough->m_descriptor.chIndex;
228  DetectorRegionIndex region = road.seed->hough->m_descriptor.region;
229  ATH_MSG_DEBUG("analyse: Seeding new road: eta maxima "
230  << road.maxima.size() << " phi " << road.phiMaxima.size() << " seed : sector " << sector << " "
231  << regionName(region) << " " << layerName(layer)
232  << " maximum " << road.seed->max << " position " << road.seed->pos << " angle " << road.seed->theta);
233 
234  if (road.phiMaxima.empty())
235  unassociatedEtaMaxima.push_back(road.maxima);
236  else {
237  for (auto& max : road.mergedPhiMaxima) { phiEtaAssMap[&max] = road.maxima; }
238  }
239  createPatternCombinations(phiEtaAssMap, *patternCombis);
240  createPatternCombinations(unassociatedEtaMaxima, *patternCombis);
241  }
242 
243  } else {
244  // now that the full hough transform is filled, order sectors by maxima
245  std::vector<HoughDataPerSector*> sectorData(state.houghDataPerSectorVec->vec.size());
246  for (unsigned int i = 0; i < state.houghDataPerSectorVec->vec.size(); ++i) sectorData[i] = &state.houghDataPerSectorVec->vec[i];
247  std::stable_sort(sectorData.begin(), sectorData.end(), SortHoughDataPerSector());
248 
249  std::vector<HoughDataPerSector*>::iterator spit = sectorData.begin();
250  std::vector<HoughDataPerSector*>::iterator spit_end = sectorData.end();
251  for (; spit != spit_end; ++spit) {
252  // get data for this sector
253  HoughDataPerSector& houghData = **spit;
254  using namespace MuonStationIndex;
255 
256  // loop over regions
257  for (int reg = 0; reg < toInt(DetectorRegionIndex::DetectorRegionIndexMax); ++reg) {
258  DetectorRegionIndex region = static_cast<DetectorRegionIndex>(reg);
259 
260  // only run analysis on sectors with maxima
261  if (houghData.nlayersWithMaxima[reg] == 0) continue;
262  ATH_MSG_DEBUG("Analyzing sector "
263  << (*spit)->sector << " " << regionName(region) << " nmax " << (*spit)->maxEtaHits()
264  << " layers with eta maxima " << houghData.nlayersWithMaxima[toInt(region)] << " hits "
265  << houghData.nmaxHitsInRegion[toInt(region)] << " layers with phi maxima "
266  << houghData.nphilayersWithMaxima[toInt(region)] << " hits " << houghData.nphimaxHitsInRegion[toInt(region)]);
267 
268  // look for maxima in the overlap regions of sectors
269  associateMaximaInNeighbouringSectors(houghData, state.houghDataPerSectorVec->vec);
270 
271  // layers in this region
272  constexpr int nlayers = toInt(LayerIndex::LayerIndexMax);
273 
274  // first link phi maxima with eta maxima
275  RegionMaximumVec unassociatedEtaMaxima(nlayers);
276  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec> phiEtaAssociations;
277  associateMaximaToPhiMaxima(region, houghData, phiEtaAssociations, unassociatedEtaMaxima);
278 
279  // create pattern combinations for combined patterns
280  createPatternCombinations(phiEtaAssociations, *patternCombis);
281 
282  // create pattern combinations for unassociated patterns
283  createPatternCombinations(unassociatedEtaMaxima, *patternCombis);
284  }
285  }
286  }
287 
288  ATH_MSG_DEBUG("Found " << patternCombis->size() << " pattern combinations " << std::endl << m_printer->print(*patternCombis));
289 
290  if (msgLvl(MSG::DEBUG)) {
291  ATH_MSG_DEBUG("Hough performance ");
292  printTruthSummary(state.truthHits, state.foundTruthHits);
293  ATH_MSG_DEBUG("Association performance ");
294  printTruthSummary(state.foundTruthHits, state.outputTruthHits);
295  }
296 
297  return {std::move(patternCombis), std::move(state.houghDataPerSectorVec)};
298  }
299 
300  void MuonLayerHoughTool::buildRoads(MaximumVec& seedMaxima, MuonHough::MuonDetectorHough& detectorHoughTransforms,
301  std::unique_ptr<HoughDataPerSectorVec>& houghDataPerSectorVec,
302  std::vector<MuonLayerHoughTool::Road>& roads) const {
303  // sort maxima according to hits
304  std::stable_sort(seedMaxima.begin(), seedMaxima.end(),
305  [](const std::shared_ptr<MuonHough::MuonLayerHough::Maximum>& m1,
306  const std::shared_ptr<MuonHough::MuonLayerHough::Maximum>& m2) { return m1->max > m2->max; });
307  // loop over seed maxima (which are maxima) that pass certain thresholds detailed in cut_values
308  std::set<std::shared_ptr<MuonHough::MuonLayerHough::Maximum>> associatedMaxima;
309  for (const auto& seed : seedMaxima) {
310  // if this maximum is already in the set of associated maxima, do not do anything
311  if (associatedMaxima.count(seed)) continue;
312 
313  // maximum becomes our new seed
314 
315  // decomposing the locality information for the seed
316  int sector = seed->hough->m_descriptor.sector;
317  MuonStationIndex::ChIndex chIndex = seed->hough->m_descriptor.chIndex;
319  using DetRegIdx = DetRegIdx;
320  DetRegIdx region = seed->hough->m_descriptor.region;
321 
322  // creating new road with said seed
323  Road road(seed);
324  ATH_MSG_DEBUG(" New seed: sector " << seed->hough->m_descriptor.sector << " " << MuonStationIndex::regionName(region)
325  << " " << MuonStationIndex::layerName(layer) << " maximum " << seed->max
326  << " position " << seed->pos << " angle " << seed->theta << " ptr " << seed.get());
331  MuonHough::HitVec::const_iterator ref_itr = std::find_if(
332  seed->hits.begin(), seed->hits.end(), [](const std::shared_ptr<MuonHough::Hit>& hit) -> bool { return hit->prd; });
333 
334  const bool isNSW = ref_itr != seed->hits.end() &&
335  (m_idHelperSvc->issTgc((*ref_itr)->prd->identify()) || m_idHelperSvc->isMM((*ref_itr)->prd->identify()));
336  // extend seed within the current sector
337  // sector indices have an offset of -1 because the numbering of the sectors are from 1 to 16 but the indices in the vertices are
338  // of course 0 to 15
339  extendSeed(detectorHoughTransforms, road, houghDataPerSectorVec->vec[sector - 1]);
340 
341  // look for maxima in the overlap regions of sectors
342  int sectorN = sector - 1;
343  if (sectorN < 1) sectorN = 16;
344  int sectorP = sector + 1;
345  if (sectorP > 16) sectorP = 1;
346 
347  // associate the road with phi maxima
348  associatePhiMaxima(road, houghDataPerSectorVec->vec[sector - 1].phiMaxVec[toInt(region)]);
349  //
350  if (m_addSectors && isNSW) {
351  extendSeed(detectorHoughTransforms, road, houghDataPerSectorVec->vec[sectorN - 1]);
352  associatePhiMaxima(road, houghDataPerSectorVec->vec[sectorN - 1].phiMaxVec[toInt(region)]);
353  extendSeed(detectorHoughTransforms, road, houghDataPerSectorVec->vec[sectorP - 1]);
354  associatePhiMaxima(road, houghDataPerSectorVec->vec[sectorP - 1].phiMaxVec[toInt(region)]);
355  }
356 
357  if (road.neighbouringRegion != DetRegIdx::DetectorRegionUnknown) {
358  associatePhiMaxima(road, houghDataPerSectorVec->vec[sector - 1].phiMaxVec[toInt(road.neighbouringRegion)]);
359  }
360  // if close to a sector boundary, try adding maxima in that sector as well
361  if (road.neighbouringSector != -1) {
362  ATH_MSG_DEBUG(" Adding neighbouring sector " << road.neighbouringSector);
363  extendSeed(detectorHoughTransforms, road,
364  houghDataPerSectorVec->vec[road.neighbouringSector - 1]);
365  associatePhiMaxima(road, houghDataPerSectorVec->vec[road.neighbouringSector - 1].phiMaxVec[toInt(region)]);
366  }
367 
368  // finally deal with the case that we have both neighbouring region and sector
369  if (road.neighbouringRegion != DetRegIdx::DetectorRegionUnknown && road.neighbouringSector != -1) {
370  associatePhiMaxima(road, houghDataPerSectorVec->vec[road.neighbouringSector - 1].phiMaxVec[toInt(road.neighbouringRegion)]);
371  }
372 
373  // merge phi maxima
374  mergePhiMaxima(road);
375 
376  // add maxima to seed exclusion list
377  associatedMaxima.insert(road.maxima.begin(), road.maxima.end());
378 
379  if (msgLevel(MSG::DEBUG)) {
380  ATH_MSG_DEBUG(" New road " << road.maxima.size());
381  for (const auto& max : road.maxima) {
382  MuonStationIndex::ChIndex chIndex = max->hough->m_descriptor.chIndex;
384  DetRegIdx region = max->hough->m_descriptor.region;
385  ATH_MSG_DEBUG(" Sector " << max->hough->m_descriptor.sector << " " << MuonStationIndex::regionName(region) << " "
386  << MuonStationIndex::layerName(layer) << " maximum " << max->max << " position "
387  << max->pos << " angle " << max->theta << " ptr " << max);
388  }
389  }
390  bool insert = true;
391  for (auto& oldRoad : roads) {
392  std::vector<std::shared_ptr<MuonHough::MuonLayerHough::Maximum>> intersection;
393  std::set_intersection(oldRoad.maximumSet.begin(), oldRoad.maximumSet.end(), road.maximumSet.begin(), road.maximumSet.end(),
394  std::back_inserter(intersection));
395  unsigned int intersectionSize = intersection.size();
396  unsigned int oldRoadSize = oldRoad.maximumSet.size();
397  unsigned int roadSize = road.maximumSet.size();
398  ATH_MSG_VERBOSE(" Overlap check " << intersectionSize << " old " << oldRoadSize << " new " << roadSize << " old ptr "
399  << oldRoad.seed);
400  if (intersectionSize == 0) continue;
401  if (intersectionSize == roadSize) {
402  insert = false; // discard
403  break;
404  } else if (intersectionSize == oldRoadSize) {
405  oldRoad = road; // replace
406  insert = false;
407  break;
408  }
409  }
410 
411  // add road to list
412  if (insert) roads.push_back(road);
413  }
414  }
415 
417  // input -> list of phiMaxima on road
418  // returns some mergedPhiMaxima -> is this "summed" over layers?
419 
420  auto maximaSortingLambda = [road](const std::shared_ptr<MuonHough::MuonPhiLayerHough::Maximum>& m1,
421  const std::shared_ptr<MuonHough::MuonPhiLayerHough::Maximum>& m2) {
422  if (m1->max != m2->max) return m1->max > m2->max;
423  // prefer the same sector as the seed sector
424  if (m1->sector != m2->sector) return m1->sector == road.seed->hough->m_descriptor.sector;
425 
426  if (m1->hits.size() != m2->hits.size()) return m1->hits.size() < m2->hits.size(); // least hits -> most collimated maximum
427 
428  if (m1->pos != m2->pos) return m1->pos < m2->pos;
429 
430  if (std::abs(m1->binposmax - m1->binposmin) == std::abs(m2->binposmax - m2->binposmin)) {
431  return (m1->binposmin) < (m2->binposmin);
432  }
433  return std::abs(m1->binposmax - m1->binposmin) < std::abs(m2->binposmax - m2->binposmin);
434  };
435 
436  std::stable_sort(road.phiMaxima.begin(), road.phiMaxima.end(), maximaSortingLambda);
437 
438  ATH_MSG_VERBOSE("Merging phi maxima " << road.phiMaxima.size());
439  std::set<MuonHough::MuonPhiLayerHough::Maximum*> associatedPhiMaxima;
440  for (auto pit = road.phiMaxima.begin(); pit != road.phiMaxima.end(); ++pit) { // loop over phi maxima
441  if (associatedPhiMaxima.count((*pit).get())) continue; // check if maximum is already in associatedPhiMaxima
442  associatedPhiMaxima.insert((*pit).get());
443  MuonHough::MuonPhiLayerHough::Maximum phiMaximum = **pit;
444  ATH_MSG_VERBOSE(" phi maxima " << phiMaximum.pos << " val " << phiMaximum.max);
445 
446  bool wasExtended = false;
447  for (auto pit1 = pit + 1; pit1 != road.phiMaxima.end(); ++pit1) {
448  if ((*pit1)->binposmax >= phiMaximum.binposmin && (*pit1)->binposmin <= phiMaximum.binposmax) {
449  ATH_MSG_VERBOSE(" merging maxima " << phiMaximum.pos << " val " << phiMaximum.max << " " << (*pit1)->pos << " val "
450  << (*pit1)->max);
451  phiMaximum.hits.insert(phiMaximum.hits.end(), (*pit1)->hits.begin(), (*pit1)->hits.end());
452  associatedPhiMaxima.insert((*pit1).get());
453  wasExtended = true;
454  }
455  }
456 
457  if (wasExtended) {
458  // refind maximum
459  MuonHough::MuonPhiLayerHough localHough(
460  60, -M_PI, M_PI, ((*pit)->hough ? (*pit)->hough->m_region : DetRegIdx::DetectorRegionUnknown));
461  MuonHough::PhiHitVec hits = phiMaximum.hits;
462  /* too ambiguous producing irreproducibilities because of sorting by pointer value
463  std::stable_sort(hits.begin(),hits.end(),[]( const MuonHough::PhiHit* h1,
464  const MuonHough::PhiHit* h2 ){ return h1->layer < h2->layer; } );
465  */
466 
467  std::stable_sort(hits.begin(), hits.end(),
468  [](const std::shared_ptr<MuonHough::PhiHit>& h1, const std::shared_ptr<MuonHough::PhiHit>& h2) {
469  if (h1->layer != h2->layer) return h1->layer < h2->layer;
470  if (h1->w != h2->w) return h1->w > h2->w;
471  if (h1->r != h2->r) return h1->r < h2->r;
472 
473  const double dPhi1 = std::abs(h1->phimax - h1->phimin);
474  const double dPhi2 = std::abs(h2->phimax - h2->phimin);
475  if (dPhi1 != dPhi2) return dPhi1 < dPhi2;
476  if (h1->phimin == h2->phimin) return h1->phimax < h2->phimax;
477  return h1->phimin < h2->phimin;
478  });
479 
480  ATH_MSG_VERBOSE(" updating phi maximum " << phiMaximum.pos << " bin " << phiMaximum.binpos << " val " << phiMaximum.max
481  << " number of hits " << hits.size());
482  if (msgLvl(MSG::VERBOSE)) localHough.setDebug(true);
483  localHough.fillLayer2(hits);
484  localHough.findMaximum(phiMaximum, 0.9);
485  localHough.associateHitsToMaximum(phiMaximum, hits);
486  ATH_MSG_VERBOSE(" updated phi maxima " << phiMaximum.pos << " bin " << phiMaximum.binpos << " val " << phiMaximum.max
487  << " number of hits " << phiMaximum.hits.size());
488  phiMaximum.hough = (*pit)->hough; // set back pointer to transform
489  }
490  road.mergedPhiMaxima.push_back(phiMaximum);
491  }
492  }
493 
494  // maximum in middle layer
495  // says look in other layers
496  // if yes, combine them
497  // gets on road
498  // roads are combinations of maxima
499 
502  MuonLayerHoughTool::HoughDataPerSector& sectorData) const { // const {
503  if (!road.seed) return;
504 
505  RegionMaximumVec& maxVec = sectorData.maxVec;
506 
507  // gather locality information on seed
509  MuonStationIndex::LayerIndex seedLayer = MuonStationIndex::toLayerIndex(seed.hough->m_descriptor.chIndex);
510  DetRegIdx region = seed.hough->m_descriptor.region;
511 
512  // loop over layers in the same region as the seed ( inner, middle, outer)
513  for (int lay = 0; lay < toInt(LayerIndex::LayerIndexMax); ++lay) {
515  if (layer == seedLayer && seed.hough->m_descriptor.sector == sectorData.sector) continue;
516 
517  // untrue -> look in neighboring layer
518  // true -> look only in this layer
519  double distanceCut = layer == seedLayer ? 500. : (double)m_extrapolationDistance;
520 
521  unsigned int layerHash = MuonStationIndex::sectorLayerHash(region, layer);
522 
523  // fetching vector of maxima for given region and layer
524  const MaximumVec& maxima = maxVec[layerHash];
525  if (maxima.empty()) continue;
526 
527  ATH_MSG_DEBUG("Associating maxima in " << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
528  << " size " << maxima.size());
529  // loop over maxima in layer
530  for (const auto& candMaximum : maxima) {
531  // extrapolate seed to layer assuming a pointing straight line or parabolic
532  // add maximum to road if close enough
533  float yloc_diff = MuonHough::extrapolate(seed, *candMaximum, m_doParabolicExtrapolation);
534  if (std::abs(MuonHough::extrapolate(seed, *candMaximum, m_doParabolicExtrapolation)) < distanceCut) {
535  ATH_MSG_VERBOSE(" Adding maximum position " << candMaximum->pos << " intersect diff" << yloc_diff);
536  road.add(candMaximum);
537  } else {
538  ATH_MSG_VERBOSE(" Maximum position: y "
539  << candMaximum->pos << " x " << candMaximum->hough->m_descriptor.referencePosition << " seed y "
540  << seed.hough->m_descriptor.referencePosition << " x " << seed.pos << " intersect diff " << yloc_diff);
541  }
542  }
543  }
544 
545  // check if the maximum is close to the detector boundary, if yes look for maxima in the neighbouring region, skip BarrelExtended
546  if (seedLayer == LayerIndex::BarrelExtended) return;
547 
548  ATH_MSG_DEBUG("Checking Barrel/Endcap overlaps: min dist edge "
549  << seed.pos - seed.hough->m_descriptor.yMinRange << " max dist edge " << seed.pos - seed.hough->m_descriptor.yMaxRange
550  << " pos " << seed.pos << " range " << seed.hough->m_descriptor.yMinRange << " "
551  << seed.hough->m_descriptor.yMaxRange);
552 
553  if (std::abs(seed.pos - seed.hough->m_descriptor.yMinRange) < 4000. ||
554  std::abs(seed.pos - seed.hough->m_descriptor.yMaxRange) < 4000.) {
555  // asumes region is barrel and looks in adjacent regions (clever logic TM here)
556  DetRegIdx neighbourRegion = DetRegIdx::Barrel;
557  if (region == DetRegIdx::Barrel) {
558  if (seed.pos < 0)
559  neighbourRegion = DetRegIdx::EndcapC;
560  else
561  neighbourRegion = DetRegIdx::EndcapA;
562  } // in all other cases the neigbourRegion is definitely barrel
563 
564  // looping over all layers in neigbouring region
565  for (int lay = 0; lay < toInt(LayerIndex::LayerIndexMax); ++lay) {
567 
568  // skip barrel combinations with BEE
569  if (region == DetRegIdx::Barrel && layer == LayerIndex::BarrelExtended) continue;
570 
571  double distanceCut = 1000.;
572 
573  // get maxima from neigboring region
574  unsigned int layerHash = MuonStationIndex::sectorLayerHash(neighbourRegion, layer);
575  const MaximumVec& maxima = maxVec[layerHash];
576  if (maxima.empty()) continue;
577  ATH_MSG_DEBUG("Associating maxima in neighbouring region " << MuonStationIndex::regionName(neighbourRegion) << " "
578  << MuonStationIndex::layerName(layer) << " hash " << layerHash
579  << " size " << maxima.size());
580 
581  // loop over maxima per layer
582  for (const auto& candMaximum : maxima) {
583  // extrapolate seed to layer assuming a pointing straight line, swap coordinates
584  float yloc_diff = MuonHough::extrapolate(seed, *candMaximum, m_doParabolicExtrapolation);
585  ATH_MSG_VERBOSE(" Maximum position: y "
586  << candMaximum->pos << " x " << candMaximum->hough->m_descriptor.referencePosition << " seed y "
587  << seed.hough->m_descriptor.referencePosition << " x " << seed.pos << " intersect diff " << yloc_diff);
588 
589  if (std::abs(yloc_diff) < distanceCut) {
590  road.add(candMaximum);
591  road.neighbouringRegion = neighbourRegion;
592  }
593  }
594  }
595  }
596 
597  // search for phiMaxima using the etaMaximum of the road in the current sector
598  std::set<const TgcClusterObj3D*> tgcClusters;
599  std::set<Identifier> triggerLayers;
600  const MaximumVec& maxima = road.maxima;
601  for (const auto& maximum : maxima) {
602  if (maximum->hough->m_descriptor.sector != sectorData.sector)
603  continue; // skip cases where a maximum on the road does not belong to the currently examined sector
604 
605  // gather tgcClusters associated to the hits of the maxima
606  for (auto ehit = maximum->hits.begin(); ehit != maximum->hits.end(); ++ehit) {
607  const MuonHough::Hit& etaHit = **ehit;
608  if (etaHit.tgc) {
609  if (!etaHit.tgc->phiCluster.empty()) tgcClusters.insert(etaHit.tgc);
610  } else if (etaHit.prd) {
611  triggerLayers.insert(m_idHelperSvc->gasGapId(etaHit.prd->identify()));
612  }
613  }
614  }
615 
616  MuonHough::MuonPhiLayerHough& phiHough =
617  detectorHoughTransforms.phiHough(region); // get phi transform in the same region as the seed
618 
619  // gather phiHits in sector that match the etaHits of the maximum
620  PhiHitVec phiHitsInMaximum;
621  PhiHitVec& phiHits = sectorData.phiHitVec[toInt(region)];
622  for (const auto& phiHit : phiHits) {
623  if (phiHit->tgc) {
624  if (tgcClusters.find(phiHit->tgc) != tgcClusters.end()) phiHitsInMaximum.push_back(phiHit);
625  } else if (phiHit->prd) {
626  if (triggerLayers.find(m_idHelperSvc->gasGapId(phiHit->prd->identify())) != triggerLayers.end())
627  phiHitsInMaximum.push_back(phiHit);
628  }
629  }
630 
631  // fill phi hits
632  ATH_MSG_DEBUG("extendSeed: Filling s" << sectorData.sector << " " << MuonStationIndex::regionName(region) << " phiHitsInMaxima "
633  << phiHitsInMaximum.size() << " phi hits: " << phiHits.size());
634 
635  if (!findMaxima(phiHough, phiHitsInMaximum, sectorData.phiMaxVec[toInt(region)], sectorData.sector) ||
636  sectorData.phiMaxVec[toInt(region)].empty()) {
637  ATH_MSG_DEBUG("extendSeed: No phi maxima found in s" << sectorData.sector << " " << MuonStationIndex::regionName(region));
638  return;
639  }
640 
641  ++sectorData.nphilayersWithMaxima[toInt(region)];
642  sectorData.nphimaxHitsInRegion[toInt(region)] += sectorData.phiMaxVec[toInt(region)].front()->max;
643 
644  ATH_MSG_DEBUG("extendSeed: Sector phiMaxima Summary: s" << sectorData.sector << " " << MuonStationIndex::regionName(region) << " "
645  << sectorData.nphilayersWithMaxima[toInt(region)]
646  << " -> nPhiMaxima: " << sectorData.phiMaxVec[toInt(region)].size()
647  << " max sum: " << sectorData.nphimaxHitsInRegion[toInt(region)]);
648  }
649 
650  // phi hits are not separated into inner middle outer
651  // maxima found in road
653  ATH_MSG_DEBUG("associateMaximaToPhiMaxima: phi maxima " << phiMaxima.size());
654  if (!road.seed) return;
655 
656  // loop over phi maxima
657  for (const auto& pmaximum : phiMaxima) {
658  // reference to phi maximum
659 
660  ATH_MSG_DEBUG(" new phi maximum " << pmaximum->max << " hits " << pmaximum->hits.size());
661 
662  // precalculate the layers + TGC clusters and add them to a set for easy access
663  std::map<Identifier, std::pair<float, float>> triggerLayersPhiMinMax;
664  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>> tgcClusters;
665 
666  // loop over hits
667  for (const auto& phiHit : pmaximum->hits) {
668  // two cases
669  // case 1: phiHit measured in TGC -> get phiHits from phiCluster
670  // case 2: phiHit is prepared raw data -> use phiHit to extend the triggerLayersPhinMinMax map
671  if (phiHit->tgc) {
672  if (phiHit->tgc->phiCluster.empty())
673  ATH_MSG_WARNING(" TGC 3D cluster without phi hits ");
674  else
675  tgcClusters[m_idHelperSvc->stationIndex(phiHit->tgc->phiCluster.front()->identify())].insert(phiHit->tgc);
676  } else if (phiHit->prd) {
677  Identifier gpId = m_idHelperSvc->gasGapId(phiHit->prd->identify());
678  auto mit = triggerLayersPhiMinMax.find(gpId);
679  if (mit == triggerLayersPhiMinMax.end())
680  triggerLayersPhiMinMax[gpId] = std::make_pair(phiHit->phimin, phiHit->phimax);
681  else {
682  mit->second.first = std::min(phiHit->phimin, mit->second.first);
683  mit->second.second = std::max(phiHit->phimax, mit->second.second);
684  }
685  }
686  }
687  // print out information on the triggerLayersPhiMinMax
688  if (msgLevel(MSG::VERBOSE)) {
689  ATH_MSG_DEBUG("Trigger layers " << triggerLayersPhiMinMax.size() << " tgc layers " << tgcClusters.size());
690  for (auto tgcit = triggerLayersPhiMinMax.begin(); tgcit != triggerLayersPhiMinMax.end(); ++tgcit) {
691  ATH_MSG_VERBOSE(" " << m_idHelperSvc->toString(tgcit->first));
692  }
693 
694  // loop over the stations and the contained tgcClusters found in the previous step, print out information
695  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>>::const_iterator stit = tgcClusters.begin();
696  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>>::const_iterator stit_end = tgcClusters.end();
697  for (; stit != stit_end; ++stit) {
698  std::set<const TgcClusterObj3D*>::const_iterator ttit = stit->second.begin();
699  std::set<const TgcClusterObj3D*>::const_iterator ttit_end = stit->second.end();
700  for (; ttit != ttit_end; ++ttit) {
701  ATH_MSG_VERBOSE(" " << m_idHelperSvc->toString((*ttit)->phiCluster.front()->identify()) << " nhits "
702  << (*ttit)->phiCluster.size());
703  }
704  }
705  }
706 
707  // check if there are matching maxima in neighbouring sectors, add maximum values if confirmation is found
708  // overlap counters
709  int noverlaps = 0;
710  int nNoOverlaps = 0;
711  float phimin{10}, phimax{-10};
712 
713  // loop over all maxima found on road
714  for (const auto& road_max : road.maxima) {
715  // get station information for maximum on road
716  MuonStationIndex::StIndex stIndex = MuonStationIndex::toStationIndex(road_max->hough->m_descriptor.chIndex);
717 
718  // loop over eta hits
719  for (const auto& etaHit : road_max->hits) {
720  if (etaHit->tgc) {
721  if (etaHit->tgc->etaCluster.empty())
722  ATH_MSG_WARNING(" TGC 3D cluster without eta hits ");
723  else {
724  if (tgcClusters[stIndex].count(etaHit->tgc)) {
725  // now loop over phi maximum and find phi hit
726  for (const auto& phiHit : pmaximum->hits) {
727  if (phiHit->tgc == etaHit->tgc) {
728  phimin = std::min(phiHit->phimin, phimin);
729  phimax = std::max(phiHit->phimax, phimax);
730  break;
731  }
732  }
733  ++noverlaps;
734  } else {
735  ++nNoOverlaps;
736  }
737  }
738  } else if (etaHit->prd) {
739  if (!m_idHelperSvc->isRpc(etaHit->prd->identify()) && !m_idHelperSvc->issTgc(etaHit->prd->identify())) continue;
740  Identifier gpId = m_idHelperSvc->gasGapId(etaHit->prd->identify());
741  auto mit = triggerLayersPhiMinMax.find(gpId);
742  if (mit == triggerLayersPhiMinMax.end())
743  ++nNoOverlaps;
744  else {
745  phimin = std::min(mit->second.first, phimin);
746  phimax = std::max(mit->second.second, phimax);
747  ++noverlaps;
748  }
749  }
750  } // loop over hits in maximum
751  } // loop over maxima in road
752 
753  // if overlaps are found, add the phi maximum in question to the road
754  if (noverlaps > 0) {
755  road.add(pmaximum);
756  // check if we are close to a sector boundary
757  std::vector<int> sectors;
758  m_sectorMapping.getSectors(phimin, sectors);
759  if (sectors.size() > 1) {
760  for (const int& sec : sectors) {
761  if (sec != road.seed->hough->m_descriptor.sector) road.neighbouringSector = sec;
762  }
763  } else {
764  std::vector<int> sectors;
765  m_sectorMapping.getSectors(phimax, sectors);
766  if (sectors.size() > 1) {
767  for (const int& sec : sectors) {
768  if (sec != road.seed->hough->m_descriptor.sector) road.neighbouringSector = sec;
769  }
770  }
771  }
772  }
773  ATH_MSG_DEBUG(" Overlap with Phi maximum: overlap " << noverlaps << " no overlap " << nNoOverlaps << " phimin " << phimin
774  << " phimax " << phimax << " neighbouring sector "
775  << road.neighbouringSector);
776  }
777  }
778 
779  // takes the maxima from a given sector and tries to associate it with the maxima of the adjacent sectors
782  std::vector<MuonLayerHoughTool::HoughDataPerSector>& houghDataPerSectorVec) const {
783  ATH_MSG_DEBUG(" looping over eta maxima");
784 
785  // now loop over eta maxima per layer
786  for (unsigned int regLay = 0; regLay < houghData.maxVec.size(); ++regLay) {
787  MaximumVec& maxima = houghData.maxVec[regLay];
788  int sector = houghData.sector;
789 
790  // loop over two neighbouring sectors
791  for (int i = 0; i < 2; ++i) {
792  // calculate neighbouring sector index
793  int sectorN = (i == 0) ? sector - 1 : sector + 1;
794  if (i == 0 && sector == 1) sectorN = 16;
795  if (i == 1 && sector == 16) sectorN = 1;
796 
797  MuonLayerHoughTool::HoughDataPerSector& houghDataN = houghDataPerSectorVec[sectorN - 1];
798 
799  MaximumVec& maximaN = houghDataN.maxVec[regLay];
800 
801  // loop over maxima in layer
802  for (const auto& maximum : maxima) {
803  // reference to maximum
804 
805  if (!maximum->hough) {
806  ATH_MSG_WARNING("Maximum without associated hough transform! ");
807  continue;
808  }
809 
810  // loop over maxima per layer in neighbouring sector
811  for (const auto& maximumN : maximaN) {
812  // reference to maximum
813  if (!maximumN->hough) {
814  ATH_MSG_WARNING("Maximum without associated hough transform! ");
815  continue;
816  }
817 
818  // calculate the position of the first maximum in the reference frame of the other sector
819  double rcor = maximumN->hough->m_descriptor.referencePosition *
820  m_sectorMapping.transformRToNeighboringSector(maximum->pos, sector, sectorN) /
821  maximum->hough->m_descriptor.referencePosition;
822  double dist = rcor - maximumN->pos;
823  ATH_MSG_DEBUG(" maximumN->hough " << maximumN->hough->m_descriptor.referencePosition << " maximum->hough "
824  << maximum->hough->m_descriptor.referencePosition << " maximumN->pos "
825  << maximumN->pos << " maximum->pos " << maximum->pos << rcor << " distance "
826  << dist);
827  if (std::abs(dist) > 100) continue;
828  houghData.maxAssociationMap[maximum.get()].push_back(maximumN);
829  houghDataN.associatedToOtherSector.insert(maximumN.get());
830 
831  ATH_MSG_DEBUG(" Found maximum in neighbouring sector: max " << maximum->max << " pos " << rcor << " maxN "
832  << maximumN->max << " pos " << maximumN->pos
833  << " distance " << dist);
834 
835  // loop over first and second maximum
836  for (int nn = 0; nn < 2; ++nn) {
837  // setting info for the debug-info objects of the hits
838  const auto& maxi = nn == 0 ? maximum : maximumN;
839  const auto& maxi2 = nn == 0 ? maximumN : maximum;
840  ATH_MSG_VERBOSE(" Maximum " << nn << " hits " << maxi->hits.size());
841  for (auto& hit : maxi->hits) {
842  if (hit->debugInfo()) {
843  hit->debugInfo()->phn = maxi2->max;
844  Identifier id = hit->tgc ? hit->tgc->etaCluster.front()->identify() : hit->prd->identify();
845  ATH_MSG_VERBOSE(" " << m_idHelperSvc->toString(id) << " setphn " << hit->debugInfo()->phn);
846  }
847  }
848  }
849  }
850  }
851  }
852  }
853  }
854 
857  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MuonLayerHoughTool::MaximumVec>& phiEtaAssociations,
858  MuonLayerHoughTool::RegionMaximumVec& unassEtaMaxima) const {
859  std::set<std::shared_ptr<MuonHough::MuonLayerHough::Maximum>> associatedMaxima;
860 
861  PhiMaximumVec& phiMaxima = houghData.phiMaxVec[MuonStationIndex::toInt(region)];
862 
863  ATH_MSG_DEBUG("associateMaximaToPhiMaxima in sector " << houghData.sector << ": phi maxima " << phiMaxima.size()); // !!!!
864  // loop over phi maxima
865  for (const auto& phiMaximum : phiMaxima) {
866  // reference to phi maximum
867 
868  ATH_MSG_DEBUG(" Considering phi maximum " << phiMaximum->max << " hits " << phiMaximum->hits.size()); // !!!!
869 
870  // store associated maxima
871  MaximumVec associatedMaximaVec; // !!!!
872 
873  // precalculate the layers + TGC clusters and add them to a set for easy access
874  // std::map< Identifier,std::pair<double,double> > triggerLayersP;
875  std::set<Identifier> triggerLayers;
876  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>> tgcClusters;
877 
878  // loop over hits
879  for (const auto& phiHit : phiMaximum->hits) {
880  if (phiHit->tgc) {
881  if (phiHit->tgc->phiCluster.empty())
882  ATH_MSG_WARNING(" TGC 3D cluster without phi hits ");
883  else
884  tgcClusters[m_idHelperSvc->stationIndex(phiHit->tgc->phiCluster.front()->identify())].insert(phiHit->tgc);
885  } else if (phiHit->prd) {
886  Identifier colId = phiHit->prd->identify();
887  Identifier layId = m_idHelperSvc->gasGapId(colId);
888  triggerLayers.insert(layId);
889  }
890  }
891  if (msgLvl(MSG::DEBUG)) {
892  ATH_MSG_DEBUG("Trigger layers " << triggerLayers.size() << " tgc layers " << tgcClusters.size());
893  for (const Identifier& id : triggerLayers) { ATH_MSG_VERBOSE(" " << m_idHelperSvc->toString(id)); }
894 
895  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>>::const_iterator stit = tgcClusters.begin();
896  std::map<MuonStationIndex::StIndex, std::set<const TgcClusterObj3D*>>::const_iterator stit_end = tgcClusters.end();
897  for (; stit != stit_end; ++stit) {
898  std::set<const TgcClusterObj3D*>::const_iterator ttit = stit->second.begin();
899  std::set<const TgcClusterObj3D*>::const_iterator ttit_end = stit->second.end();
900  for (; ttit != ttit_end; ++ttit) {
901  ATH_MSG_VERBOSE(" " << m_idHelperSvc->toString((*ttit)->phiCluster.front()->identify()) << " nhits "
902  << (*ttit)->phiCluster.size());
903  }
904  }
905  }
906 
907  ATH_MSG_DEBUG(" looping over eta maxima");
908 
909  // now loop over eta maxima per layer
910  for (unsigned int lay = 0; lay < toInt(LayerIndex::LayerIndexMax); ++lay) {
912  unsigned int layerHash = MuonStationIndex::sectorLayerHash(region, layer);
913  MaximumVec& maxima = houghData.maxVec[layerHash];
914  if (maxima.empty()) continue;
916 
917  // loop over maxima per layer
918  for (const auto& maximum : maxima) {
919  // skip maxima that were already associated to a neighbouring sector
920  if (houghData.associatedToOtherSector.count(maximum.get())) continue;
921 
922  // check if there are matching maxima in neighbouring sectors, add maximum values if confirmation is found
923  float totmax = 0;
924  int ntrigconfirm = 0;
925  MaximumAssociationMap::iterator pos = houghData.maxAssociationMap.find(maximum.get());
926  if (pos != houghData.maxAssociationMap.end()) {
927  for (const auto& max_itr : pos->second) {
928  totmax = std::max(max_itr->max, totmax);
929  ntrigconfirm += max_itr->triggerConfirmed;
930  }
931  }
932  totmax += maximum->max;
933  ntrigconfirm += maximum->triggerConfirmed;
934 
935  ATH_MSG_DEBUG(" new eta maximum " << MuonStationIndex::stName(stIndex) << " val " << maximum->max
936  << " neightbour confirmed value " << totmax << " trigger confirmations "
937  << ntrigconfirm);
938 
939  // overlap counters
940  int nmmHits{0}, ntgcOverlaps{0}, nrpcOverlaps{0}, nstgcOverlaps{0}, nstgcNoOverlaps{0};
941 
942  // loop over hits
943  for (const auto& etaHit : maximum->hits) {
944  if (etaHit->tgc) {
945  if (tgcClusters[stIndex].count(etaHit->tgc))
946  ++ntgcOverlaps;
947 
948  } else if (etaHit->prd) {
949  Identifier layId = m_idHelperSvc->gasGapId(etaHit->prd->identify());
950  ATH_MSG_VERBOSE(" eta layer hit " << m_idHelperSvc->toString(layId));
951  if (m_idHelperSvc->isMM(layId)) ++nmmHits;
952  if (triggerLayers.count(layId)) {
953  if (m_idHelperSvc->isRpc(layId))
954  ++nrpcOverlaps;
955  else if (m_idHelperSvc->issTgc(layId))
956  ++nstgcOverlaps;
957  } else {
958  if (m_idHelperSvc->issTgc(layId))
959  ++nstgcNoOverlaps;
960  }
961  }
962  }
963 
964  // cuts on NSW endcap only for now
965  if (nmmHits + nstgcNoOverlaps + nstgcOverlaps > 0) {
966  // select
967  if (maximum->pos < 1200.) {
968  if (totmax < 8) {
969  ATH_MSG_DEBUG(" maximum failed cut " << totmax << " cut 8, position " << maximum->pos);
970  continue;
971  }
972  } else if (maximum->pos > 4300.) {
973  if (totmax < 8) {
974  ATH_MSG_DEBUG(" maximum failed cut " << totmax << " cut 8, position " << maximum->pos);
975  continue;
976  }
977  } else {
978  if (totmax < 12) {
979  ATH_MSG_DEBUG(" maximum failed cut " << totmax << " cut 12, position " << maximum->pos);
980  continue;
981  }
982  }
983  }
984 
985  ATH_MSG_DEBUG(" Overlap with Phi maximum: tgc " << ntgcOverlaps << " stgc " << nstgcOverlaps << " rpc " << nrpcOverlaps
986  << " nphiTgc " << tgcClusters[stIndex].size() << " trigLay "
987  << triggerLayers.size());
988  if (stIndex == StIndex::EM && !tgcClusters[stIndex].empty() && ntgcOverlaps == 0) {
989  ATH_MSG_VERBOSE(" No association in StationLayer " << MuonStationIndex::stName(stIndex) << " tgcs overlaps "
990  << ntgcOverlaps << " on phi maximum "
991  << tgcClusters[stIndex].size());
992  continue;
993  }
994  if (stIndex == StIndex::EI && !tgcClusters[stIndex].empty() && ntgcOverlaps == 0) {
995  ATH_MSG_VERBOSE(" No association in StationLayer " << MuonStationIndex::stName(stIndex) << " tgcs overlaps "
996  << ntgcOverlaps << " on phi maximum "
997  << tgcClusters[stIndex].size());
998  continue;
999  }
1000  if (stIndex == StIndex::EI && nstgcOverlaps == 0 && nstgcNoOverlaps != 0) {
1001  ATH_MSG_VERBOSE(" No association in StationLayer " << MuonStationIndex::stName(stIndex)
1002  << " stgcs without overlaps " << nstgcNoOverlaps);
1003  continue;
1004  }
1005  // require STGC confirmation
1006  if (m_requireTriggerConfirmationNSW && nmmHits > 0 && ntrigconfirm == 0) continue;
1007 
1008  associatedMaxima.insert(maximum);
1009  associatedMaximaVec.push_back(maximum);
1010 
1011  // check if there are matching maxima in neighbouring sectors
1012  if (pos != houghData.maxAssociationMap.end()) {
1013  associatedMaxima.insert(pos->second.begin(), pos->second.end());
1014  associatedMaximaVec.insert(associatedMaximaVec.end(), pos->second.begin(), pos->second.end());
1015  }
1016  }
1017  }
1018 
1019  if (associatedMaximaVec.empty()) continue;
1020  ATH_MSG_DEBUG(" processed phi maximum, associated eta maxima " << associatedMaximaVec.size());
1021  phiEtaAssociations[phiMaximum.get()] = std::move(associatedMaximaVec);
1022  }
1023 
1024  // finally idenitify all unassociated maxima and add them to the unassociated maxima list
1025  // now loop over eta maxima per layer
1026  for (unsigned int lay = 0; lay < toInt(LayerIndex::LayerIndexMax); ++lay) {
1028  unsigned int layerHash = MuonStationIndex::sectorLayerHash(region, layer);
1029 
1030  if (lay >= unassEtaMaxima.size()) {
1031  ATH_MSG_WARNING(" size of unassEtaMaxima too small for region " << unassEtaMaxima.size() << " region "
1032  << MuonStationIndex::regionName(region));
1033  break;
1034  }
1035  MaximumVec& maxima = houghData.maxVec[layerHash];
1036 
1037  // loop over maxima per layer
1038  for (const auto& mit : maxima) {
1039  if (associatedMaxima.count(mit)) continue;
1040  unassEtaMaxima[lay].push_back(mit);
1041  ATH_MSG_DEBUG(" unassociated maximum in layer " << MuonStationIndex::layerName(layer) << " max-val " << mit->max);
1042  }
1043  }
1044  }
1045 
1047  MuonPatternCombinationCollection& patternCombis) const {
1048  ATH_MSG_DEBUG("Creating pattern combinations for eta patterns ");
1049 
1050  std::vector<MuonPatternChamberIntersect> chamberData;
1051 
1052  // bool isEndcap = maxima.size() == 5;
1053 
1054  // loop over layers
1055  for (const auto& max_sec : maxima) {
1056  // create vector for prds per chamber
1057  std::map<Identifier, std::set<const Trk::PrepRawData*>> prdsPerChamber;
1058 
1059  // loop over maxima per layer
1060  for (const auto& max : max_sec) {
1061  ATH_MSG_DEBUG(" new maximum " << max->max << " hits " << max->hits.size());
1062 
1063  // sanity check
1064  if (max->hits.empty()) {
1065  ATH_MSG_WARNING(" Maximum without hits ");
1066  continue;
1067  }
1068  ATH_MSG_DEBUG(" adding hits " << max->hits.size());
1069 
1070  // loop over hits in maximum and add them to the hit list
1071  for (const auto& hit : max->hits) {
1072  if (hit->tgc) {
1073  const Identifier chId = m_idHelperSvc->chamberId(hit->tgc->etaCluster.front()->identify());
1074  prdsPerChamber[chId].insert(hit->tgc->etaCluster.begin(), hit->tgc->etaCluster.end());
1075  } else if (hit->prd) {
1076  const Identifier chId = m_idHelperSvc->chamberId(hit->prd->identify());
1077  prdsPerChamber[chId].insert(hit->prd);
1078  }
1079  }
1080  }
1081 
1082  auto sortPrdIds = [](const Trk::PrepRawData* prd1, const Trk::PrepRawData* prd2) {
1083  return prd1->identify() < prd2->identify();
1084  };
1085  std::map<Identifier, std::set<const Trk::PrepRawData*>>::iterator chit = prdsPerChamber.begin();
1086  std::map<Identifier, std::set<const Trk::PrepRawData*>>::iterator chit_end = prdsPerChamber.end();
1087  for (; chit != chit_end; ++chit) {
1088  ATH_MSG_DEBUG("Adding chamber " << m_idHelperSvc->toStringChamber(chit->first) << " hits " << chit->second.size());
1089  std::vector<const Trk::PrepRawData*> prds;
1090  prds.insert(prds.end(), chit->second.begin(), chit->second.end());
1091  std::stable_sort(prds.begin(), prds.end(), sortPrdIds);
1092  const Trk::PrepRawData& prd = **prds.begin();
1093  Amg::Vector3D gpos = prd.detectorElement()->surface(prd.identify()).center();
1094  // create intersection and add it to combination
1095  ATH_MSG_DEBUG("Adding chamber with intersect phi direction " << gpos.phi() << " theta " << gpos.theta());
1096  MuonPatternChamberIntersect intersect(gpos, gpos.unit(), prds);
1097  chamberData.push_back(intersect);
1098  }
1099  }
1100  if (chamberData.empty()) return;
1101 
1102  MuonPatternCombination* combi = new MuonPatternCombination(nullptr, chamberData);
1103 
1104  ATH_MSG_DEBUG(" creating new unassociated " << m_printer->print(*combi));
1105  patternCombis.push_back(combi);
1106  }
1107 
1109  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MuonLayerHoughTool::MaximumVec>& phiEtaAssociations,
1110  MuonPatternCombinationCollection& patternCombis) const {
1111  ATH_MSG_DEBUG("Creating pattern combinations from eta/phi combinations " << phiEtaAssociations.size());
1112 
1113  // loop over the phi maxima
1114  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec>::const_iterator pit = phiEtaAssociations.begin();
1115  std::map<MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec>::const_iterator pit_end = phiEtaAssociations.end();
1116  for (; pit != pit_end; ++pit) {
1117  if (pit->second.empty()) continue;
1118 
1119  // collect phi hits per chamber
1120  std::map<Identifier, std::set<const Trk::PrepRawData*>> phiHitsPerChamber;
1121 
1122  // loop over hits
1123  for (const auto& hit : pit->first->hits) {
1124  if (hit->tgc) {
1125  const Identifier chId = m_idHelperSvc->chamberId(hit->tgc->phiCluster.front()->identify());
1126  phiHitsPerChamber[chId].insert(hit->tgc->phiCluster.begin(), hit->tgc->phiCluster.end());
1127  } else if (hit->prd) {
1128  const Identifier chId = m_idHelperSvc->chamberId(hit->prd->identify());
1129  phiHitsPerChamber[chId].insert(hit->prd);
1130  }
1131  }
1132 
1133  // create chamber intersections
1134  std::vector<MuonPatternChamberIntersect> chamberData;
1135  std::set<Identifier> addedPhiHits;
1136 
1137  // create vector for prds per chamber
1138  std::map<Identifier, std::set<const Trk::PrepRawData*>> prdsPerChamber;
1139 
1140  // store position and direction of the first maximum in the chamber layer
1141  std::map<MuonStationIndex::ChIndex, std::pair<Amg::Vector3D, Amg::Vector3D>> directionsPerChamberLayer;
1142 
1143  // loop over eta maxima
1144  for (const auto& max : pit->second) {
1145  ATH_MSG_DEBUG(" new maximum " << max->max << " hits " << max->hits.size());
1146 
1147  if (!max->hough) { ATH_MSG_WARNING("Maximum without associated Hough Transform"); }
1148 
1149  // sanity check
1150  if (max->hits.empty()) {
1151  ATH_MSG_WARNING(" Maximum without hits ");
1152  continue;
1153  }
1154  ATH_MSG_DEBUG(" adding hits " << max->hits.size());
1155 
1156  // loop over hits in maximum and add them to the hit list
1157  for (const auto& hit : max->hits) {
1158  Identifier chId;
1159  if (hit->tgc) {
1160  chId = m_idHelperSvc->chamberId(hit->tgc->etaCluster.front()->identify());
1161  prdsPerChamber[chId].insert(hit->tgc->etaCluster.begin(), hit->tgc->etaCluster.end());
1162  } else if (hit->prd) {
1163  chId = m_idHelperSvc->chamberId(hit->prd->identify());
1164  prdsPerChamber[chId].insert(hit->prd);
1165  } else {
1166  ATH_MSG_WARNING("Hit without associated PRDs");
1167  continue;
1168  }
1169  // the first time we have a maximun in this layer store the position and direction
1170  MuonStationIndex::ChIndex chIndex = m_idHelperSvc->chamberIndex(chId);
1171  if (!directionsPerChamberLayer.count(chIndex)) {
1172  // eta maximum has z(r) and theta parameters but these are local
1173  double maxpos = max->pos;
1174  double refPlane = 0.;
1175  bool isBarrel = !m_idHelperSvc->isEndcap(chId) || chIndex == ChIndex::BEE;
1176  if (max->hough)
1177  refPlane = max->hough->m_descriptor.referencePosition;
1178  else if (hit->tgc)
1179  refPlane = hit->tgc->getEdge(TgcEdge::LowEtaLowPhi).z();
1180  else if (isBarrel)
1181  refPlane = hit->prd->detectorElement()->surface(hit->prd->identify()).center().perp();
1182  else
1183  refPlane = hit->prd->detectorElement()->surface(hit->prd->identify()).center().z();
1184 
1185  double r = isBarrel ? refPlane : maxpos;
1186  double z = isBarrel ? maxpos : refPlane;
1187  double theta = max->theta;
1188 
1189  // go to global
1190  double sign = 1.;
1191  if (isBarrel) {
1192  theta += M_PI_2;
1193  sign = -1.;
1194  }
1195 
1196  // phi maximum has one phi from position assume global Phi definition
1197  double phi = pit->first->pos; // phiCor(pit->first->pos,pit->first->sector,false);
1198 
1199  CxxUtils::sincos scphi(phi);
1200  double sinphi = scphi.sn;
1201  double cosphi = scphi.cs;
1202 
1203  CxxUtils::sincos sctheta(theta);
1204  double sintheta = sctheta.sn;
1205  double costheta = sctheta.cs;
1206 
1207  std::pair<Amg::Vector3D, Amg::Vector3D>& posDir = directionsPerChamberLayer[chIndex];
1208  posDir.first = Amg::Vector3D(r * cosphi, r * sinphi, z);
1209  posDir.second = Amg::Vector3D(sign * cosphi * costheta, sign * sinphi * costheta, sintheta);
1211  << " setting position: perp " << posDir.first.perp() << " z " << posDir.first.z() << " phi pos "
1212  << posDir.first.phi() << " direction phi " << posDir.second.phi() << " theta pos "
1213  << posDir.first.theta() << " direction theta " << posDir.second.theta() << " ref perp " << r << " z "
1214  << z << " phi " << phi << " theta " << theta);
1215  if (posDir.first.dot(posDir.second) < 0.) {
1216  ATH_MSG_WARNING(" direction not pointing to IP " << posDir.first.unit().dot(posDir.second));
1217  }
1218  }
1219 
1220  std::map<Identifier, std::set<const Trk::PrepRawData*>>::iterator pos = phiHitsPerChamber.find(chId);
1221  if (pos != phiHitsPerChamber.end()) {
1222  std::pair<std::set<Identifier>::iterator, bool> ipos = addedPhiHits.insert(chId);
1223  if (ipos.second) { prdsPerChamber[chId].insert(pos->second.begin(), pos->second.end()); }
1224  }
1225  }
1226  }
1227 
1228  auto sortPrdIds = [](const Trk::PrepRawData* prd1, const Trk::PrepRawData* prd2) {
1229  return prd1->identify() < prd2->identify();
1230  };
1231  std::map<Identifier, std::set<const Trk::PrepRawData*>>::iterator chit = prdsPerChamber.begin();
1232  std::map<Identifier, std::set<const Trk::PrepRawData*>>::iterator chit_end = prdsPerChamber.end();
1233  for (; chit != chit_end; ++chit) {
1234  ATH_MSG_DEBUG("Adding chamber " << m_idHelperSvc->toStringChamber(chit->first) << " hits " << chit->second.size());
1235  std::vector<const Trk::PrepRawData*> prds;
1236  prds.insert(prds.end(), chit->second.begin(), chit->second.end());
1237  std::stable_sort(prds.begin(), prds.end(), sortPrdIds);
1238  const Trk::PrepRawData& prd = **prds.begin();
1239 
1240  MuonStationIndex::ChIndex chIndex = m_idHelperSvc->chamberIndex(prd.identify());
1241  std::map<MuonStationIndex::ChIndex, std::pair<Amg::Vector3D, Amg::Vector3D>>::const_iterator pos =
1242  directionsPerChamberLayer.find(chIndex);
1245  if (pos != directionsPerChamberLayer.end()) {
1246  gpos = pos->second.first;
1247  gdir = pos->second.second;
1248  } else {
1249  ATH_MSG_WARNING("No global position and direction found, calculating from surface");
1250  gpos = prd.detectorElement()->surface(prd.identify()).center();
1251  gdir = -1 * gpos.unit();
1252  }
1253 
1254  ATH_MSG_DEBUG("Creating intersection " << MuonStationIndex::chName(chIndex) << " setting position: perp " << gpos.perp()
1255  << " z " << gpos.z() << " phi pos " << gpos.phi() << " direction phi " << gdir.phi()
1256  << " theta pos " << gpos.theta() << " theta " << gdir.theta() << " hits "
1257  << prds.size());
1258 
1259  // create intersection and add it to combination
1260  MuonPatternChamberIntersect intersect(gpos, gdir, prds);
1261  chamberData.push_back(intersect);
1262 
1263  }
1264  if (chamberData.empty()) continue;
1265  if (addedPhiHits.empty()) {
1266  ATH_MSG_DEBUG("No phi hits selected, skipping combi ");
1267  continue;
1268  }
1269  MuonPatternCombination* combi = new MuonPatternCombination(nullptr, chamberData);
1270  ATH_MSG_DEBUG("adding pattern combination with chambers " << chamberData.size() << " phi layers " << addedPhiHits.size()
1271  << std::endl
1272  << m_printer->print(*combi));
1273  patternCombis.push_back(combi);
1274  }
1275  }
1276 
1279  MuonLayerHoughTool::MaximumVec& maxima) const {
1280  if (hits.empty()) return false;
1281 
1282  using namespace MuonStationIndex;
1283  if (hough.m_descriptor.chIndex == ChIndex::ChUnknown || hough.m_descriptor.chIndex == ChIndex::ChIndexMax) {
1284  Identifier id = hits.front()->tgc ? hits.front()->tgc->etaCluster.front()->identify() : hits.front()->prd->identify();
1285  ATH_MSG_WARNING("Bad ChIndex " << m_idHelperSvc->toString(id) << " " << chName(hough.m_descriptor.chIndex));
1286  return false;
1287  }
1288 
1289  // populate hough transform with hits
1290  std::stable_sort(hits.begin(), hits.end(), MuonHough::SortHitsPerLayer());
1291  if (m_debugHough) hough.setDebug(true);
1292  hough.fillLayer2(hits);
1293 
1294  Identifier id_hit = hits.front()->tgc ? hits.front()->tgc->etaCluster.front()->identify() : hits.front()->prd->identify();
1295  MuonHough::MuonLayerHoughSelector selectorLoose;
1297 
1298  if (m_idHelperSvc->issTgc(id_hit) || m_idHelperSvc->isMM(id_hit)) {
1299  selectorLoose = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 3.9)});
1300  selector = MuonHough::MuonLayerHoughSelector({std::make_pair(0, 7.9)});
1301  } else {
1302  selectorLoose = m_selectorsLoose[toInt(hough.m_descriptor.chIndex)];
1303  selector = m_selectors[toInt(hough.m_descriptor.chIndex)];
1304  }
1305 
1306  // MuonStationIndex::StIndex stIndex = MuonStationIndex::toStationIndex(hough.m_descriptor.chIndex);
1307  unsigned int nmaxima = 0;
1308  while (nmaxima < 5) {
1310  if (hough.findMaximum(maximum, selectorLoose)) {
1311  hough.associateHitsToMaximum(maximum, hits);
1312  ATH_MSG_VERBOSE("findMaxima: Found Eta Maximum "
1313  << nmaxima << " " << maximum.max << " trigConfirmed " << maximum.triggerConfirmed << " pos " << maximum.pos
1314  << " theta " << maximum.theta << " binPos " << maximum.binpos << " binRange " << maximum.binposmin << " -- "
1315  << maximum.binposmax << " binTheta " << maximum.bintheta << " nHits " << maximum.hits.size());
1316 
1317  int nmdt = 0;
1318  int nmm = 0;
1319  int nstgc = 0;
1320 
1321  const unsigned int nHitsInMaximum = maximum.hits.size();
1322  for (unsigned int i = 0; i < nHitsInMaximum; ++i) {
1323  MuonHough::Hit& hit = *(maximum.hits[i]);
1324  Identifier id = hit.tgc ? hit.tgc->etaCluster.front()->identify() : hit.prd->identify();
1325  int nhits = hit.tgc ? hit.tgc->etaCluster.size() : 1;
1326 
1327  nmdt += m_idHelperSvc->isMdt(id);
1328  nstgc += m_idHelperSvc->issTgc(id);
1329  nmm += m_idHelperSvc->isMM(id);
1330 
1331  ATH_MSG_VERBOSE("findMaxima: hit " << hit.layer << " " << m_idHelperSvc->toString(id) << " hits " << nhits);
1332  }
1333 
1334  // only store maxima that have MDT hits
1335  if (nmdt > 0 || (nmm + nstgc) > 0) {
1336  maxima.emplace_back(std::make_unique<MuonHough::MuonLayerHough::Maximum>(maximum));
1337  // add to seed list if
1338  if (maximum.max > selector.getCutValue(maximum.pos)) seedMaxima.push_back(maxima.back());
1339  ++nmaxima;
1340  }
1341  hough.fillLayer2(maximum.hits, true);
1342  } else {
1343  if (nmaxima > 0) { ATH_MSG_VERBOSE("findMaxima: No more maxima found " << nmaxima); }
1344  // ?!? if nmaximo == 0 here the function should return false, I think
1345  break;
1346  }
1347  }
1348  return true;
1349  }
1350 
1352  MuonLayerHoughTool::PhiMaximumVec& maxima, int sector) const {
1353  if (hits.empty()) return false;
1354 
1355  std::stable_sort(hits.begin(), hits.end(), MuonHough::SortHitsPerLayer());
1356  if (m_debugHough) hough.setDebug(true);
1357  hough.fillLayer2(hits);
1358 
1359  unsigned int nmaxima = 0;
1360  while (nmaxima < 5) {
1362  if (hough.findMaximum(maximum, 1.9)) {
1363  hough.associateHitsToMaximum(maximum, hits);
1364 
1365  ATH_MSG_DEBUG("findMaxima(Phi): Found Phi maximum " << nmaxima << " height " << maximum.max << " pos " << maximum.pos
1366  << " bin pos " << maximum.binpos << " binRange " << maximum.binposmin
1367  << " -- " << maximum.binposmax << " nHits " << maximum.hits.size());
1368 
1369  const unsigned int nHitsInMaximum = maximum.hits.size();
1370  for (unsigned int i = 0; i < nHitsInMaximum; ++i) {
1371  MuonHough::PhiHit& hit = *(maximum.hits[i]);
1372  Identifier id = hit.tgc ? hit.tgc->phiCluster.front()->identify() : hit.prd->identify();
1373 
1374  int nhits = hit.tgc ? hit.tgc->phiCluster.size() : 1;
1375  ATH_MSG_VERBOSE("findMaxima(Phi) phiHit " << m_idHelperSvc->toString(id) << " hits " << nhits);
1376  }
1377 
1378  maximum.sector = sector; // very fragile passing on of sector
1379 
1380  // check if the maximum is already filled, if so, don't add it again
1381  bool maximum_matched = false;
1382  for (auto pit = maxima.begin(); pit != maxima.end(); ++pit) {
1383  // reference to phi maximum
1384  MuonHough::MuonPhiLayerHough::Maximum& pmaximum = **pit;
1385  if (pmaximum.sector == maximum.sector && pmaximum.max == maximum.max && pmaximum.pos == maximum.pos &&
1386  pmaximum.hits.size() == maximum.hits.size() && pmaximum.binpos == maximum.binpos &&
1387  pmaximum.binposmin == maximum.binposmin && pmaximum.binposmax == maximum.binposmax) {
1388  ATH_MSG_DEBUG("extendSeed: sector has already been added! Skip. ");
1389  bool maximum_hitmatched = true; // check if there is a hit that is not the same
1390  for (unsigned int k = 0; k < maximum.hits.size(); ++k) {
1391  if (maximum.hits[k] != pmaximum.hits[k]) { // directly compare pointer address
1392  maximum_hitmatched = false;
1393  break;
1394  }
1395  }
1396  if (maximum_hitmatched) {
1397  maximum_matched = true;
1398  break;
1399  }
1400  }
1401  }
1402  // remove the hits from hough
1403  hough.fillLayer2(maximum.hits, true);
1404  if (maximum_matched) {
1405  //++nmaxima;
1406  continue;
1407  } else {
1408  maxima.push_back(std::make_shared<MuonHough::MuonPhiLayerHough::Maximum>(maximum));
1409  ++nmaxima;
1410  }
1411  } else {
1412  if (nmaxima > 0) { ATH_MSG_VERBOSE("findMaxima(Phi): No more maxima found " << nmaxima); }
1413  // ?!? same here, the function should return false if nothing was found, right?
1414  break;
1415  }
1416  }
1417  hough.reset();
1418  return true;
1419  }
1420 
1421  void MuonLayerHoughTool::fillHitsPerSector(const EventContext& ctx, State& state, const int sector,
1422  const CollectionsPerSector& collectionsPerSector, const MdtPrepDataContainer* mdtCont,
1423  const CscPrepDataContainer* /*cscCont*/, const TgcPrepDataContainer* tgcCont,
1424  const RpcPrepDataContainer* rpcCont, const sTgcPrepDataContainer* stgcCont,
1425  const MMPrepDataContainer* mmCont) const {
1426  MuonLayerHoughTool::HoughDataPerSector& houghData = state.houghDataPerSectorVec->vec[sector - 1];
1427  houghData.sector = sector;
1428  // loop over all possible station layers in the sector
1429  for (unsigned int tech = 0; tech < m_ntechnologies; ++tech) {
1430  for (unsigned int layerHash = 0; layerHash < MuonStationIndex::sectorLayerHashMax(); ++layerHash) {
1431  const HashVec& hashes = collectionsPerSector.technologyRegionHashVecs[tech][layerHash];
1432  if (hashes.empty()) continue;
1433  auto regionLayer = MuonStationIndex::decomposeSectorLayerHash(layerHash);
1434 
1435  for (const IdentifierHash& id_hash : hashes) {
1436  // !?! else if made by Felix
1437  if (mdtCont && !mdtCont->empty() && tech == toInt(TechnologyIndex::MDT)) {
1438  const MdtPrepDataCollection* pos = mdtCont->indexFindPtr(id_hash);
1439  if (pos) fill(ctx, state.truthHits, *pos, houghData.hitVec[layerHash]);
1440  } else if (rpcCont && !rpcCont->empty() && tech == toInt(TechnologyIndex::RPC)) {
1441  const RpcPrepDataCollection* pos = rpcCont->indexFindPtr(id_hash);
1442  if (pos) fill(ctx, state.truthHits, *pos, houghData.hitVec[layerHash], houghData.phiHitVec[toInt(regionLayer.first)]);
1443  } else if (tgcCont && !tgcCont->empty() && tech == toInt(TechnologyIndex::TGC)) {
1444  const TgcPrepDataCollection* pos = tgcCont->indexFindPtr(id_hash);
1445  if (pos)
1446  fill(ctx, state.truthHits, state.houghDataPerSectorVec->tgcClusteringObjs, *pos, houghData.hitVec[layerHash],
1447  houghData.phiHitVec[toInt(regionLayer.first)], collectionsPerSector.sector);
1448  } else if (stgcCont && !stgcCont->empty() && tech == toInt(TechnologyIndex::STGC)) {
1449  const sTgcPrepDataCollection* pos = stgcCont->indexFindPtr(id_hash);
1450  if (pos)
1451  fill(ctx, state.truthHits, *pos, houghData.hitVec[layerHash], houghData.phiHitVec[toInt(regionLayer.first)],
1452  collectionsPerSector.sector);
1453  } else if (mmCont && !mmCont->empty() && tech == toInt(TechnologyIndex::MM)) {
1454  const MMPrepDataCollection* pos = mmCont->indexFindPtr(id_hash);
1455  if (pos) fill(ctx, state.truthHits, *pos, houghData.hitVec[layerHash]);
1456  }
1457  }
1458  }
1459  }
1460  }
1461 
1462  void MuonLayerHoughTool::matchTruth(std::set<Identifier>& truthHits, const PRD_MultiTruthCollection& truthCol, const Identifier& id,
1463  MuonHough::HitDebugInfo& debug) const {
1464  typedef PRD_MultiTruthCollection::const_iterator iprdt;
1465  std::pair<iprdt, iprdt> range = truthCol.equal_range(id);
1466  // Loop over particles contributing to this cluster
1467  for (iprdt i = range.first; i != range.second; ++i) {
1468  if (!i->second.isValid()) {
1469  ATH_MSG_WARNING("Unexpected invalid HepMcParticleLink in PRD_MultiTruthCollection");
1470  } else {
1471  const HepMcParticleLink& link = i->second;
1472  if (link.cptr() && abs(link.cptr()->pdg_id()) == 13) {
1473  debug.uniqueID = HepMC::uniqueID(link);
1474  debug.pdgId = link.cptr()->pdg_id();
1475  truthHits.insert(id);
1476  }
1477  }
1478  }
1479  }
1480 
1481  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits, const MdtPrepDataCollection& mdts,
1483  if (mdts.empty()) return;
1484  auto truthCollections = m_truthNames.makeHandles(ctx);
1485  Identifier chid = mdts.identify();
1486  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1487  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1488  int sector = m_idHelperSvc->sector(chid);
1489  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1490  bool barrelLike = (region == DetRegIdx::Barrel || layer == LayerIndex::BarrelExtended);
1491  unsigned int nmdts(mdts.size()), nmdtsBad{0};
1492  for (const MdtPrepData* prd : mdts) {
1493  if (prd->adc() < 50 || prd->status() != MdtStatusDriftTime) {
1494  ++nmdtsBad;
1495  continue;
1496  }
1497  const Identifier id = prd->identify();
1498  float r = rCor(*prd);
1499  float x = barrelLike ? r : prd->globalPosition().z();
1500  float y = barrelLike ? prd->globalPosition().z() : r;
1501  int sublayer = sublay(id);
1502 
1503  float ymin = y - prd->localPosition()[Trk::locR];
1504  float ymax = y + prd->localPosition()[Trk::locR];
1505  MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology, sector, region, layer, sublayer);
1506  debug->time = prd->tdc();
1507  debug->r = prd->localPosition()[Trk::locR];
1508 
1509  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1510  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1511  MuonHough::Hit* hit = new MuonHough::Hit(sublayer, x, ymin, ymax, 1., debug, prd);
1512  hits.emplace_back(hit);
1513  }
1514 
1515  ATH_MSG_DEBUG("fillMDT: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1517  << " -> hits: " << nmdts << " bad " << nmdtsBad << " isSmallChamber "
1518  << m_idHelperSvc->isSmallChamber(chid));
1519  }
1520 
1521  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits, const CscPrepDataCollection& cscs, HitVec& hits,
1522  PhiHitVec& phiHits) const {
1524  return;
1525  auto truthCollections = m_truthNames.makeHandles(ctx);
1526  Identifier chid = cscs.identify();
1527  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1528  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1529  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1530  int sector = m_idHelperSvc->sector(chid);
1531  unsigned int neta{0}, nphi{0};
1532  for (const CscPrepData* prd : cscs) {
1533  const bool meas_phi = m_idHelperSvc->rpcIdHelper().measuresPhi(prd->identify());
1534  nphi += meas_phi;
1535  neta += !meas_phi;
1536  }
1537  ATH_MSG_DEBUG("fillCscs: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1539  << " -> eta hits " << neta << " phi hits " << nphi);
1540  for (const CscPrepData* prd : cscs) {
1541  const Identifier id = prd->identify();
1542  int sublayer = sublay(id);
1543  MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology, sector, region, layer, sublayer);
1544  debug->isEtaPhi = (neta && nphi);
1545  debug->trigConfirm = 1;
1546  debug->time = prd->time();
1547  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1548  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1549  float weight = (neta && nphi) ? 2 : 1;
1550  if (m_idHelperSvc->rpcIdHelper().measuresPhi(id)) {
1551  const float r = rCor(*prd);
1552  const float phi = prd->globalPosition().phi();
1553  const double phi1 = phi; // phiCor(phi,sector);
1554  debug->r = -99999;
1555  MuonHough::PhiHit* hit = new MuonHough::PhiHit(sublayer, r, phi1, phi1, weight, debug, prd);
1556  phiHits.emplace_back(hit);
1557  } else {
1558  const float x = rCor(*prd);
1559  const float y = prd->globalPosition().z();
1560  const float stripCor = 0.5 * prd->detectorElement()->StripWidth(false);
1561  const float ymin = y - stripCor;
1562  const float ymax = y + stripCor;
1563  debug->r = stripCor;
1564  MuonHough::Hit* hit = new MuonHough::Hit(sublayer, x, ymin, ymax, weight, debug, prd);
1565  hits.emplace_back(hit);
1566  }
1567  }
1568  }
1569  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits, const RpcPrepDataCollection& rpcs,
1571  if (rpcs.empty()) return;
1572  auto truthCollections = m_truthNames.makeHandles(ctx);
1573  Identifier chid = rpcs.identify();
1574  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1575  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1576  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1577  int sector = m_idHelperSvc->sector(chid);
1578  // check whether there are eta and phi hits
1579  unsigned int neta{0}, nphi{0};
1580  for (const RpcPrepData* prd : rpcs) {
1581  const bool meas_phi = m_idHelperSvc->rpcIdHelper().measuresPhi(prd->identify());
1582  nphi += meas_phi;
1583  neta += !meas_phi;
1584  }
1585  ATH_MSG_DEBUG("fillRPC: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1587  << " -> eta hits " << neta << " phi hits " << nphi);
1588 
1589  for (const RpcPrepData* prd : rpcs) {
1590  const Identifier id = prd->identify();
1591  int sublayer = sublay(id);
1592  MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology, sector, region, layer, sublayer);
1593  debug->isEtaPhi = (neta && nphi);
1594  debug->trigConfirm = 1;
1595  debug->time = prd->time();
1596  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1597  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1598  float weight = (neta && nphi) ? 2 : 1;
1599  if (m_idHelperSvc->rpcIdHelper().measuresPhi(id)) {
1600  const float r = rCor(*prd);
1601  const float phi = prd->globalPosition().phi();
1602  const double phi1 = phi; // phiCor(phi,sector);
1603  debug->r = -99999;
1604  MuonHough::PhiHit* hit = new MuonHough::PhiHit(sublayer, r, phi1, phi1, weight, debug, prd);
1605  phiHits.emplace_back(hit);
1606  } else {
1607  const float x = rCor(*prd);
1608  const float y = prd->globalPosition().z();
1609  const float stripCor = 0.5 * prd->detectorElement()->StripWidth(false);
1610  const float ymin = y - stripCor;
1611  const float ymax = y + stripCor;
1612  debug->r = stripCor;
1613  MuonHough::Hit* hit = new MuonHough::Hit(sublayer, x, ymin, ymax, weight, debug, prd);
1614  hits.emplace_back(hit);
1615  }
1616  }
1617  }
1618 
1619  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits, const MMPrepDataCollection& mms,
1621  if (mms.empty()) return;
1622  auto truthCollections = m_truthNames.makeHandles(ctx);
1623  Identifier chid = mms.identify();
1624  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1625  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1626  int sector = m_idHelperSvc->sector(chid);
1627  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1628  ATH_MSG_DEBUG("fillMM: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1629  << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer) << " -> hits "
1630  << mms.size());
1631 
1632  std::array<double,8> multiplicity{};
1633  for (const MMPrepData* prd : mms) {
1634  const Identifier id = prd->identify();
1635  int sublayer = sublay(id) % 10;
1636  multiplicity[sublayer]++;
1637  }
1638 
1639  if( msgLvl(MSG::DEBUG) ){
1640  for (int i = 0; i<8 ; i++) if(multiplicity[i]>0) ATH_MSG_DEBUG(" sublayer " << i << " hits " << multiplicity[i]);
1641  }
1642 
1643  for (const MMPrepData* prd : mms) {
1644  const Identifier id = prd->identify();
1645  float x = prd->globalPosition().z();
1646  float y = rCor(*prd);
1647  int sublayer = sublay(id) % 10;
1648  float stripCor = prd->detectorElement()->getDesign(id)->inputPitch;
1649  float ymin = y - stripCor;
1650  float ymax = y + stripCor;
1651 
1652  //Downweight noise channels
1653  const double weight = 1. / std::max(1., multiplicity[sublayer]);
1654 
1655  MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology, sector, region, layer, sublayer);
1656  debug->r = stripCor;
1657  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1658  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1659  std::unique_ptr<MuonHough::Hit> hit = std::make_unique<MuonHough::Hit>(sublayer, x, ymin, ymax, weight, debug, prd);
1660  hits.emplace_back(std::move(hit));
1661  }
1662  }
1663 
1664  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits, const sTgcPrepDataCollection& stgcs,
1665  MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits, int selectedSector) const {
1666  if (stgcs.empty()) return;
1667  auto truthCollections = m_truthNames.makeHandles(ctx);
1668  Identifier chid = stgcs.identify();
1669  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1670  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1671  int sector = m_idHelperSvc->sector(chid);
1672  bool isNeighbouringSector = sector != selectedSector;
1673  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1674  ATH_MSG_DEBUG("fillsTGC: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1676  << " -> hits: " << stgcs.size());
1677 
1678  for (const sTgcPrepData* prd : stgcs) {
1679  const Identifier id = prd->identify();
1680  int channelType = m_idHelperSvc->stgcIdHelper().channelType(id);
1681 
1682  // only pick up phi hits in neighbouring sectors
1683  if (isNeighbouringSector && channelType == 1) continue;
1684  int sublayer = sublay(id);
1685 
1686  std::unique_ptr<MuonHough::HitDebugInfo> debug =
1687  std::make_unique<MuonHough::HitDebugInfo>(technology, sector, region, layer, sublayer);
1688  debug->isEtaPhi = 1;
1689  debug->trigConfirm = true;
1690 
1691  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1692  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1693  if (m_idHelperSvc->stgcIdHelper().channelType(id) == 1) {
1694  // eta strips
1695  float x = prd->globalPosition().z();
1696  float y = rCor(*prd);
1697  float stripCor = 1.5; // get from det el
1698  const MuonGM::MuonChannelDesign* design = prd->detectorElement()->getDesign(id);
1699  if (design) {
1700  double stripWidth = design->inputWidth;
1701  double stripLength = design->channelLength(m_idHelperSvc->stgcIdHelper().channel(id));
1702  if (m_debugHough) ATH_MSG_DEBUG(" eta strip width " << stripWidth << " stripLength " << stripLength);
1703  stripCor = 0.5 * stripWidth;
1704  }
1705  debug->r = stripCor;
1706  float ymin = y - stripCor;
1707  float ymax = y + stripCor;
1708  MuonHough::Hit* hit = new MuonHough::Hit(sublayer, x, ymin, ymax, 1., debug.release(), prd);
1709  hits.emplace_back(hit);
1710  } else {
1711  double chWidth = 0;
1713  if (m_idHelperSvc->stgcIdHelper().channelType(id) == 0) {
1714  const MuonGM::MuonPadDesign* design = prd->detectorElement()->getPadDesign(id);
1715  if (!design) {
1716  ATH_MSG_WARNING("No design found for " << m_idHelperSvc->toString(id));
1717  continue;
1718  }
1719  // get the pad width from the detector design
1720  chWidth = 0.5 * design->channelWidth(prd->localPosition(), true);
1721  ATH_MSG_DEBUG(" Pad chWidth " << chWidth << " phi global " << prd->globalPosition().phi());
1722  } else if (m_idHelperSvc->stgcIdHelper().channelType(id) == 2) {
1723  const MuonGM::MuonChannelDesign* design = prd->detectorElement()->getDesign(id);
1724  if (!design) {
1725  ATH_MSG_WARNING("No design found for " << m_idHelperSvc->toString(id));
1726  continue;
1727  }
1728  chWidth = 0.5 * design->channelWidth();
1729  ATH_MSG_DEBUG(" Wire Gang chWidth " << chWidth << " phi global " << prd->globalPosition().phi());
1730  }
1731 
1732  Amg::Vector2D lp1(prd->localPosition().x() + chWidth, prd->localPosition().y());
1733  Amg::Vector3D gp1;
1734  prd->detectorElement()->surface(id).localToGlobal(lp1, gp1, gp1);
1735 
1736  lp1[0] = prd->localPosition().x() - chWidth;
1737  Amg::Vector3D gp2;
1738  prd->detectorElement()->surface(id).localToGlobal(lp1, gp2, gp2);
1739 
1740  double phi1 = gp1.phi();
1741  double phi2 = gp2.phi();
1742  double phi1c = phi1; // phiCor(phi1,selectedSector);
1743  double phi2c = phi2; // phiCor(phi2,selectedSector);
1744  double phi_check = std::abs(xAOD::P4Helpers::deltaPhi(phi1c, phi2c));
1745  if (phi_check > 0.3) {
1746  ATH_MSG_WARNING("bad local phi: in " << phi1 << ", " << phi2 << " sector phi "
1747  << m_sectorMapping.sectorPhi(selectedSector) << " phicor " << phi1c << ", "
1748  << phi2c);
1749  }
1750  if (isNeighbouringSector &&
1751  !(m_sectorMapping.insideSector(selectedSector, phi1) || m_sectorMapping.insideSector(selectedSector, phi2))) {
1752  ATH_MSG_DEBUG("Dropping phi hit in neighbouring sector " << m_idHelperSvc->toString(id) << " phi min "
1753  << std::min(phi1c, phi2c) << " max " << std::max(phi1c, phi2c)
1754  << " global phi: in " << phi1 << ", " << phi2 << " sector phi "
1755  << m_sectorMapping.sectorPhi(selectedSector));
1756  continue;
1757  }
1758  float r = rCor(*prd);
1759  float phiMin = std::min(phi1c, phi2c);
1760  float phiMax = std::max(phi1c, phi2c);
1761  ATH_MSG_VERBOSE("Phi hit " << m_idHelperSvc->toString(id) << " r " << r << " phi min " << phiMin << " phi max "
1762  << phiMax << " bc " << debug->uniqueID << " chw " << chWidth << " trigC "
1763  << debug->trigConfirm << " g phi " << phi1 << " " << phi2);
1764  MuonHough::PhiHit* phiHit =
1765  new MuonHough::PhiHit(sublayer, r, phiMin, phiMax, 1, debug.release(), prd);
1766  phiHits.emplace_back(phiHit);
1767  }
1768  }
1769  }
1770 
1771  void MuonLayerHoughTool::fill(const EventContext& ctx, std::set<Identifier>& truthHits,
1772  std::vector<std::unique_ptr<TgcHitClusteringObj>>& tgcClusteringObjs, const TgcPrepDataCollection& tgcs,
1773  MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits, int sector) const {
1774  if (tgcs.empty()) return;
1775  tgcClusteringObjs.push_back(std::make_unique<TgcHitClusteringObj>(&m_idHelperSvc->tgcIdHelper()));
1776  TgcHitClusteringObj& clustering = *tgcClusteringObjs.back();
1777  std::vector<const TgcPrepData*> prds;
1778  prds.insert(prds.begin(), tgcs.begin(), tgcs.end());
1779  clustering.cluster(prds);
1780  clustering.buildClusters3D();
1781 
1782  Identifier chid = tgcs.identify();
1783  DetRegIdx region = m_idHelperSvc->regionIndex(chid);
1784  MuonStationIndex::LayerIndex layer = m_idHelperSvc->layerIndex(chid);
1785 
1786  if (clustering.clusters3D.empty()) {
1787  ATH_MSG_DEBUG("TgcHitClusteringObj, no 3D clusters! ");
1788  if (msgLvl(MSG::DEBUG)) {
1789  for (const TgcPrepData* prd : tgcs) { ATH_MSG_DEBUG(" " << m_idHelperSvc->toString(prd->identify())); }
1790  }
1791  return;
1792  }
1793  if (clustering.bestEtaCluster().empty()) {
1794  ATH_MSG_DEBUG("TgcHitClusteringObj, no eta cluster selected! ");
1795  if (msgLvl(MSG::DEBUG)) {
1796  for (const TgcPrepData* prd : prds) { ATH_MSG_DEBUG(" " << m_idHelperSvc->toString(prd->identify())); }
1797  }
1798  return;
1799  }
1800  auto truthCollections = m_truthNames.makeHandles(ctx);
1801  std::vector<int> sectors;
1802  getSectors(clustering.clusters3D.front(), sectors);
1803  unsigned int technology = toInt(m_idHelperSvc->technologyIndex(chid));
1804  for (unsigned int si = 0; si < sectors.size(); ++si) {
1805  if (sectors[si] != sector) continue;
1806 
1807  for (const TgcClusterObj3D& cl : clustering.clusters3D) {
1808  const Identifier id = cl.etaCluster.front()->identify();
1809 
1810  double x = cl.getEdge(TgcEdge::LowEtaLowPhi).z();
1811  double y11 = rCor(cl, TgcEdge::LowEtaLowPhi, sector);
1812  double y12 = rCor(cl, TgcEdge::LowEtaHighPhi, sector);
1813  double y21 = rCor(cl, TgcEdge::LowEtaLowPhi, sector);
1814  double y22 = rCor(cl, TgcEdge::HighEtaHighPhi, sector);
1815  double phi11 = cl.getEdge(TgcEdge::LowEtaLowPhi).phi();
1816  double phi12 = cl.getEdge(TgcEdge::LowEtaHighPhi).phi();
1817  double phi21 = cl.getEdge(TgcEdge::LowEtaLowPhi).phi();
1818  double phi22 = cl.getEdge(TgcEdge::HighEtaHighPhi).phi();
1819  double ymin = std::min(std::min(y11, y12), std::min(y21, y22));
1820  double ymax = std::max(std::max(y11, y12), std::max(y21, y22));
1821  double phimin = std::min(std::min(phi11, phi12), std::min(phi21, phi22));
1822  double phimax = std::max(std::max(phi11, phi12), std::max(phi21, phi22));
1823  double phi1 = phimin; // phiCor(phimin,sector);
1824  double phi2 = phimax; // phiCor(phimax,sector);
1825  int sublayer = sublay(id, x);
1826  ATH_MSG_VERBOSE("Cluster "<<m_idHelperSvc->toString(id)<<" x: "<<x<<", y11: "<<y11
1827  <<", y12: "<<y12<<", y21: "<<y21<<", y22: "<<y22<<", phi11: "<<phi11<<", "
1828  <<"phi12: "<<phi12<<", phi21: "<<phi21<<", phi22: "<<phi22<<" ymin: "<<ymin<<", ymax: "<<ymax
1829  <<", phimin: "<<phimin<<", phimax: "<<phimax);
1830 
1831  MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology, sector, region, layer, sublayer);
1832  debug->clusterSize = cl.etaCluster.size();
1833  debug->clusterLayers = 2;
1834  debug->isEtaPhi = true;
1835  debug->time = cl.etaCluster.front()->getBcBitMap();
1836  std::map<unsigned int, unsigned int>::const_iterator pos = m_techToTruthNameIdx.find(technology);
1837  if (pos != m_techToTruthNameIdx.end()) { matchTruth(truthHits, *truthCollections[pos->second], id, *debug); }
1838 
1840  phiDebug->clusterSize = cl.phiCluster.size();
1841  phiDebug->clusterLayers = 1;
1842  phiDebug->isEtaPhi = true;
1843 
1844  std::unique_ptr<MuonHough::Hit> hit = std::make_unique<MuonHough::Hit>(sublayer, x, ymin, ymax, 2, debug, nullptr, &cl);
1845  std::unique_ptr<MuonHough::PhiHit> phiHit =
1846  std::make_unique<MuonHough::PhiHit>(sublayer, y11, phi1, phi2, 2, phiDebug, nullptr, &cl);
1847  hits.emplace_back(std::move(hit));
1848  phiHits.emplace_back(std::move(phiHit));
1849  }
1850  }
1851  ATH_MSG_DEBUG("fillTGC: Filling " << m_idHelperSvc->toStringChamber(chid) << ": loc s" << sector << " "
1853  << " -> etaHits: " << hits.size() << " phiHits: " << phiHits.size()
1854  << " sectors: " << sectors.size());
1855  }
1856 
1858  insertHash(m_idHelperSvc->sector(id), hash, id);
1859  }
1860 
1861  void MuonLayerHoughTool::insertHash(int sector, const IdentifierHash& hash, const Identifier& id) const{
1862  MuonStationIndex::TechnologyIndex techIndex = m_idHelperSvc->technologyIndex(id);
1863  int sectorLayerHash = MuonStationIndex::sectorLayerHash(m_idHelperSvc->regionIndex(id), m_idHelperSvc->layerIndex(id));
1864  m_collectionsPerSector[sector - 1].technologyRegionHashVecs[toInt(techIndex)][sectorLayerHash].push_back(hash);
1865  }
1866 
1867  // all chambers are mapped onto a layer and sector map
1868  void MuonLayerHoughTool::initializeSectorMapping(const EventContext& ctx) const{
1869  if (m_sectorSetup) return;
1870  std::lock_guard kuchen(m_mutex);
1871  // cppcheck-suppress identicalConditionAfterEarlyExit; false positive
1872  if (m_sectorSetup) return;
1873  SG::ReadCondHandle<MuonGM::MuonDetectorManager> detMgr{m_muonManagerKey, ctx};
1874  m_collectionsPerSector.resize(MuonStationIndex::numberOfSectors());
1875  // set sector numbers
1876  unsigned int nsectorHashMax = MuonStationIndex::sectorLayerHashMax();
1877  for (unsigned int i = 0; i < m_collectionsPerSector.size(); ++i) {
1878  m_collectionsPerSector[i].sector = i + 1;
1879  m_collectionsPerSector[i].technologyRegionHashVecs.resize(m_ntechnologies);
1880  for (auto it = m_collectionsPerSector[i].technologyRegionHashVecs.begin();
1881  it != m_collectionsPerSector[i].technologyRegionHashVecs.end(); ++it) {
1882  it->resize(nsectorHashMax);
1883  }
1884  }
1885  ATH_MSG_DEBUG("Initializing hashes: number of sectors " << MuonStationIndex::numberOfSectors() << " technologies "
1886  << m_ntechnologies << " sectorLayers "
1888  // loop over all available MDT collection identifiers and order them per sector
1889 
1890  auto loadHashes = [this] (const MuonIdHelper& idHelper){
1891  auto it = idHelper.module_begin();
1892  const auto it_end = idHelper.module_end();
1893  for (; it != it_end; ++it) {
1895  idHelper.get_module_hash(*it, hash);
1896  insertHash(hash, *it);
1897  }
1898  };
1899 
1900  if (m_idHelperSvc->hasMDT()) {
1901  loadHashes(m_idHelperSvc->mdtIdHelper());
1902  }
1903  if (m_idHelperSvc->hasRPC()) {
1904  loadHashes(m_idHelperSvc->rpcIdHelper());
1905  }
1906  if (m_idHelperSvc->hasCSC()) {
1907  loadHashes(m_idHelperSvc->cscIdHelper());
1908  }
1909  // loop over all available MM collection identifiers and order them per sector
1910  if (m_idHelperSvc->hasMM()) {
1911  auto it = m_idHelperSvc->mmIdHelper().detectorElement_begin();
1912  const auto it_end = m_idHelperSvc->mmIdHelper().detectorElement_end();
1913  for (; it != it_end; ++it) {
1915  m_idHelperSvc->mmIdHelper().get_module_hash(*it, hash);
1916  insertHash(hash, *it);
1917  }
1918  }
1919  // loop over all available STGC collection identifiers and order them per sector
1920  if (m_idHelperSvc->hasSTGC()) {
1921  auto it = m_idHelperSvc->stgcIdHelper().detectorElement_begin();
1922  const auto it_end = m_idHelperSvc->stgcIdHelper().detectorElement_end();
1923  for (; it != it_end; ++it) {
1925  m_idHelperSvc->stgcIdHelper().get_module_hash(*it, hash);
1926  int sector = m_idHelperSvc->sector(*it);
1927  insertHash(sector, hash, *it);
1928  int sectorU = sector != 1 ? sector - 1 : 16;
1929  int sectorD = sector != 16 ? sector + 1 : 1;
1930  insertHash(sectorU, hash, *it);
1931  insertHash(sectorD, hash, *it);
1932  }
1933  }
1934 
1935  if (m_idHelperSvc->hasTGC()) {
1936  // loop over all available TGC collection identifiers and order them per sector
1937  auto it = m_idHelperSvc->tgcIdHelper().module_begin();
1938  const auto it_end = m_idHelperSvc->tgcIdHelper().module_end();
1939  for (; it != it_end; ++it) {
1940  const MuonGM::TgcReadoutElement* detEl = detMgr->getTgcReadoutElement(*it);
1942  m_idHelperSvc->tgcIdHelper().get_module_hash(*it, hash);
1943  int nstrips = detEl->nStrips(1);
1944  const Amg::Vector3D p1 = detEl->channelPos(1, 1, 1);
1945  const Amg::Vector3D p2 = detEl->channelPos(1, 1, nstrips);
1946  std::vector<int> sectors1{}, sectors2{};
1947  getSectors(p1, sectors1);
1948  getSectors(p2, sectors2);
1949  std::unordered_set<int> added{};
1950  for (const int sector : sectors1) {
1951  insertHash(sector, hash, *it);
1952  added.insert(sector);
1953  }
1954  for (const int sector: sectors2) {
1955  if (added.insert(sector).second){
1956  insertHash(sector, hash, *it);
1957  }
1958  }
1959  }
1960  }
1961 
1962 
1963  ATH_MSG_DEBUG(" Printing collections per sector, number of technologies " << m_ntechnologies);
1964  for (int sector = 1; sector <= 16; ++sector) {
1965  DetRegIdx currentRegion = DetRegIdx::DetectorRegionUnknown;
1966  ATH_MSG_DEBUG(" sector " << sector);
1967  TechnologyRegionHashVec& vec = m_collectionsPerSector[sector - 1].technologyRegionHashVecs;
1968  for (unsigned int hash = 0; hash < nsectorHashMax; ++hash) {
1969  std::pair<DetRegIdx, MuonStationIndex::LayerIndex> regionLayer =
1971 
1972  if (regionLayer.first != currentRegion) ATH_MSG_DEBUG(" " << MuonStationIndex::regionName(regionLayer.first));
1973  bool first = true;
1974  currentRegion = regionLayer.first;
1975  for (unsigned int tech = 0; tech < m_ntechnologies; ++tech) {
1976  std::stable_sort(vec[tech][hash].begin(), vec[tech][hash].end());
1977  if (!vec[tech][hash].empty()) {
1978  if (msgLvl(MSG::DEBUG)) {
1979  if (first) {
1980  ATH_MSG_DEBUG(" " << std::setw(7) << MuonStationIndex::layerName(regionLayer.second));
1981  first = false;
1982  }
1983  ATH_MSG_DEBUG(" " << std::setw(4)
1985  << " " << std::setw(4) << vec[tech][hash].size());
1986  }
1987  }
1988  }
1989  }
1990  }
1991  m_sectorSetup = true;
1992  }
1993 
1994  void MuonLayerHoughTool::printTruthSummary(std::set<Identifier>& truth, std::set<Identifier>& found) const {
1995  if (truth.size() == found.size()) {
1996  ATH_MSG_DEBUG(" All hits found: truth " << truth.size() << " found " << found.size());
1997  } else {
1998  ATH_MSG_DEBUG(" Some truth hits not found: truth " << truth.size() << " found " << found.size());
1999  std::vector<Identifier> result(truth.size() - found.size());
2001  std::set_difference(truth.begin(), truth.end(), found.begin(), found.end(), result.begin());
2002  result.resize(pos - result.begin());
2003  for (std::vector<Identifier>::iterator it = result.begin(); it != result.end(); ++it) {
2004  ATH_MSG_DEBUG(" " << m_idHelperSvc->toString(*it));
2005  }
2006  }
2007  }
2008 } // namespace Muon
python.root_lsr_rank.hashes
hashes
Definition: root_lsr_rank.py:34
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Muon::MuonStationIndex::LayerIndex
LayerIndex
enum to classify the different layers in the muon spectrometer
Definition: MuonStationIndex.h:38
beamspotman.r
def r
Definition: beamspotman.py:674
MuonGM::MuonPadDesign
Parameters defining the design of the readout sTGC pads.
Definition: MuonPadDesign.h:40
Muon::MuonLayerHoughTool::associateMaximaToPhiMaxima
void associateMaximaToPhiMaxima(MuonStationIndex::DetectorRegionIndex region, HoughDataPerSector &houghData, std::map< MuonHough::MuonPhiLayerHough::Maximum *, MaximumVec > &phiEtaAssociations, std::vector< MaximumVec > &unassEtaMaxima) const
Definition: MuonLayerHoughTool.cxx:855
Muon::MuonPrepDataContainer
Template for Muon PRD containers (which are basically collections of MuonPrepDataCollections).
Definition: MuonPrepDataContainer.h:42
ymin
double ymin
Definition: listroot.cxx:63
MuonHough::PhiHitVec
std::vector< std::shared_ptr< MuonHough::PhiHit > > PhiHitVec
Definition: MuonPhiLayerHough.h:20
MuonHough::MuonPhiLayerHough::Maximum::pos
float pos
Definition: MuonPhiLayerHough.h:27
Muon::MMPrepData
Class to represent MM measurements.
Definition: MMPrepData.h:22
GenEvent.h
MuonPatternChamberIntersect.h
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
get_generator_info.result
result
Definition: get_generator_info.py:21
Muon::TgcHitClusteringObj::cluster
bool cluster(const HitList &col)
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
Muon::HoughDataPerSec::phiHitVec
RegionPhiHitVec phiHitVec
Definition: HoughDataPerSec.h:53
Muon::HoughDataPerSec::phiMaxVec
RegionPhiMaximumVec phiMaxVec
Definition: HoughDataPerSec.h:55
Muon::MuonLayerHoughTool::HitVec
HoughDataPerSec::HitVec HitVec
Definition: MuonLayerHoughTool.h:59
Muon::MuonStationIndex::ChIndex::EEL
@ EEL
Muon::MuonStationIndex::sectorLayerHashMax
constexpr unsigned int sectorLayerHashMax()
maximum create a hash out of region and layer
Definition: MuonStationIndex.h:110
Muon::MuonStationIndex::ChIndex::EML
@ EML
Muon::MuonStationIndex::ChIndex::EOS
@ EOS
sTgcReadoutElement.h
Muon::MuonStationIndex::TechnologyIndex::RPC
@ RPC
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
index
Definition: index.py:1
Muon::MuonStationIndex::StIndex::EM
@ EM
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
PRD_MultiTruthCollection
A PRD is mapped onto all contributing particles.
Definition: PRD_MultiTruthCollection.h:24
MuonHough::MuonPhiLayerHough::Maximum::hough
const MuonPhiLayerHough * hough
Definition: MuonPhiLayerHough.h:37
ExaTrkXUtils::buildRoads
std::vector< std::vector< vertex_t > > buildRoads(const Graph &G, vertex_t starting_node, std::function< std::vector< vertex_t >(const Graph &, vertex_t)> next_node_fn, std::map< vertex_t, bool > &used_hits)
Definition: ExaTrkXUtils.cxx:243
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
TruthParticleContainer.h
Muon::MuonStationIndex::sectorLayerHash
unsigned int sectorLayerHash(DetectorRegionIndex detectorRegionIndex, LayerIndex layerIndex)
create a hash out of region and layer
Muon::MuonStationIndex::ChIndex::BIL
@ BIL
MuonHough::HitDebugInfo
struct containing additional debug information on the hits that is not needed for the actual alg but ...
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:26
MuonGM::MuonChannelDesign::inputWidth
double inputWidth
Definition: MuonChannelDesign.h:36
MuonHough::MuonPhiLayerHough::Maximum::binposmax
int binposmax
Definition: MuonPhiLayerHough.h:31
Muon::MuonStationIndex::stName
const std::string & stName(StIndex index)
convert StIndex into a string
Definition: MuonStationIndex.cxx:104
Muon::MuonStationIndex::LayerIndex::BarrelExtended
@ BarrelExtended
EE.
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
Muon::MuonStationIndex::TechnologyIndex
TechnologyIndex
enum to classify the different layers in the muon spectrometer
Definition: MuonStationIndex.h:54
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
Muon::MuonLayerHoughTool::RegionMaximumVec
HoughDataPerSec::RegionMaximumVec RegionMaximumVec
Definition: MuonLayerHoughTool.h:66
skel.it
it
Definition: skel.GENtoEVGEN.py:407
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Muon::MuonLayerHoughTool::Road::seed
std::shared_ptr< MuonHough::MuonLayerHough::Maximum > seed
Definition: MuonLayerHoughTool.h:78
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
MuonHough::HitDebugInfo::clusterLayers
int clusterLayers
cluster size
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:45
Muon::MuonLayerHoughTool::initialize
virtual StatusCode initialize() override
Definition: MuonLayerHoughTool.cxx:23
MuonHough::MuonLayerHough
Definition: MuonLayerHough.h:59
Muon::MuonLayerHoughTool::extendSeed
void extendSeed(MuonHough::MuonDetectorHough &detectorHoughTransforms, Road &road, HoughDataPerSector &sectorData) const
Definition: MuonLayerHoughTool.cxx:500
xAOD::P4Helpers::deltaPhi
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
Definition: xAODP4Helpers.h:69
MuonHough::MuonLayerHoughSelector
Definition: MuonLayerHoughSelector.h:12
MuonHough::MuonPhiLayerHough::fillLayer2
void fillLayer2(const PhiHitVec &hits, bool substract=false) const
Definition: MuonPhiLayerHough.cxx:29
Muon::HoughDataPerSec::nphilayersWithMaxima
std::array< int, detRegMax > nphilayersWithMaxima
Definition: HoughDataPerSec.h:57
Muon::HoughDataPerSec::nlayersWithMaxima
std::array< int, detRegMax > nlayersWithMaxima
Definition: HoughDataPerSec.h:56
Muon::HoughDataPerSec::nmaxHitsInRegion
std::array< int, detRegMax > nmaxHitsInRegion
Definition: HoughDataPerSec.h:58
vec
std::vector< size_t > vec
Definition: CombinationsGeneratorTest.cxx:9
MuonHough::MuonLayerHough::Maximum::binpos
int binpos
Definition: MuonLayerHough.h:75
Muon::MuonLayerHoughTool::State::foundTruthHits
std::set< Identifier > foundTruthHits
Definition: MuonLayerHoughTool.h:127
Trk::locR
@ locR
Definition: ParamDefs.h:44
MuonHough::HitDebugInfo::clusterSize
int clusterSize
index of reconstructed muon
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:44
Muon::MuonLayerHoughTool::findMaxima
bool findMaxima(MaximumVec &seedMaxima, MuonHough::MuonLayerHough &hough, HitVec &hits, MaximumVec &maxima) const
Definition: MuonLayerHoughTool.cxx:1277
MuonHough::Hit::layer
int layer
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:77
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Muon::MuonStationIndex::ChIndex::EIS
@ EIS
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Muon
NRpcCablingAlg reads raw condition data and writes derived condition data to the condition store.
Definition: TrackSystemController.h:45
Muon::MuonStationIndex::DetectorRegionIndex::DetectorRegionIndexMax
@ DetectorRegionIndexMax
Muon::MuonStationIndex::ChIndex::ChIndexMax
@ ChIndexMax
read_hist_ntuple.h1
h1
Definition: read_hist_ntuple.py:21
MuonHough::Hit
struct containing all hit information needed for the Hough transform
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:60
Muon::MuonLayerHoughTool::HashVec
std::vector< IdentifierHash > HashVec
Definition: MuonLayerHoughTool.h:48
x
#define x
Muon::MuonStationIndex::TechnologyIndex::MDT
@ MDT
intersection
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Definition: compareFlatTrees.cxx:25
Muon::TgcHitClusteringObj::buildClusters3D
bool buildClusters3D()
Definition: TgcHitClustering.cxx:65
Muon::MdtStatusDriftTime
@ MdtStatusDriftTime
The tube produced a vaild measurement.
Definition: MdtDriftCircleStatus.h:34
Muon::MuonStationIndex::numberOfSectors
constexpr unsigned numberOfSectors()
return total number of sectors
Definition: MuonStationIndex.h:118
MuonHough::MuonPhiLayerHough::associateHitsToMaximum
void associateHitsToMaximum(Maximum &maximum, const PhiHitVec &hits) const
Definition: MuonPhiLayerHough.cxx:246
std::stable_sort
void stable_sort(std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, std::reverse_iterator< DataModel_detail::iterator< DVL > > end, Compare comp)
Specialization of stable_sort for DataVector/List.
Definition: DVL_algorithms.h:711
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
ReadCondHandle.h
MuonHough::MuonPhiLayerHough::Maximum::binposmin
int binposmin
Definition: MuonPhiLayerHough.h:30
Muon::MuonLayerHoughTool::Road::maximumSet
std::set< std::shared_ptr< MuonHough::MuonLayerHough::Maximum > > maximumSet
Definition: MuonLayerHoughTool.h:87
Muon::MuonLayerHoughTool::TechnologyRegionHashVec
std::vector< RegionHashVec > TechnologyRegionHashVec
Definition: MuonLayerHoughTool.h:50
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
DetType::Barrel
@ Barrel
Definition: DetType.h:14
MuonGM::TgcReadoutElement::channelPos
Amg::Vector3D channelPos(const Identifier &id) const
Returns the position of the active channel (wireGang or strip)
MuonSegmentContainer.h
Muon::MuonLayerHoughTool::PhiMaximumVec
HoughDataPerSec::PhiMaximumVec PhiMaximumVec
Definition: MuonLayerHoughTool.h:64
Muon::MuonLayerHoughTool::Road::add
void add(std::shared_ptr< MuonHough::MuonLayerHough::Maximum > max)
Definition: MuonLayerHoughTool.h:79
IdentifiableContainerMT::empty
bool empty() const
return true if container is empty
Definition: IdentifiableContainerMT.h:244
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
MuonLayerHoughTool.h
Muon::MuonLayerHoughTool::analyse
std::pair< std::unique_ptr< MuonPatternCombinationCollection >, std::unique_ptr< HoughDataPerSectorVec > > analyse(State &state) const
Definition: MuonLayerHoughTool.cxx:170
Muon::HoughDataPerSec::nphimaxHitsInRegion
std::array< int, detRegMax > nphimaxHitsInRegion
Definition: HoughDataPerSec.h:59
Muon::MuonStationIndex::ChIndex::EIL
@ EIL
Muon::MuonStationIndex::toStationIndex
StIndex toStationIndex(ChIndex index)
convert ChIndex into StIndex
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
MuonHough::MuonLayerHough::Maximum::binposmin
int binposmin
Definition: MuonLayerHough.h:76
Muon::TgcHitClusteringObj::bestEtaCluster
const HitList & bestEtaCluster() const
Definition: TgcHitClustering.cxx:12
MuonHough::MuonLayerHough::Maximum
struct representing the maximum in the hough space
Definition: MuonLayerHough.h:64
Muon::MuonStationIndex::TechnologyIndex::MM
@ MM
Muon::SortHoughDataPerSector
Definition: MuonLayerHoughTool.h:218
Muon::MuonLayerHoughTool::State::seedMaxima
MaximumVec seedMaxima
Definition: MuonLayerHoughTool.h:124
MuonHough::MuonPhiLayerHough
Definition: MuonPhiLayerHough.h:22
MuonHough::MuonLayerHough::Maximum::theta
float theta
Definition: MuonLayerHough.h:69
MuonIdHelper
Definition: MuonIdHelper.h:80
Muon::HoughDataPerSec::associatedToOtherSector
std::set< MuonHough::MuonLayerHough::Maximum * > associatedToOtherSector
Definition: HoughDataPerSec.h:62
python.changerun.m1
m1
Definition: changerun.py:30
MuonHough::MuonLayerHough::Maximum::max
float max
Definition: MuonLayerHough.h:67
MuonHough::PhiHit::tgc
const Muon::TgcClusterObj3D * tgc
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:127
CxxUtils::sincos::cs
double cs
Definition: sincos.h:54
Muon::MuonStationIndex::ChIndex::BIS
@ BIS
Muon::MuonStationIndex::TechnologyIndex::TGC
@ TGC
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
Muon::RpcPrepData
Class to represent RPC measurements.
Definition: RpcPrepData.h:35
Muon::MuonPatternChamberIntersect
This class holds information needed for the Moore and MoMu pattern recognition for a muon chamber.
Definition: MuonPatternChamberIntersect.h:38
Muon::MuonStationIndex::chIndex
ChIndex chIndex(const std::string &index)
convert ChIndex name string to enum
Definition: MuonStationIndex.cxx:11
python.FPGATrackSimAnalysisConfig.hough
hough
Definition: FPGATrackSimAnalysisConfig.py:793
Muon::CscPrepData
Class representing clusters from the CSC.
Definition: CscPrepData.h:39
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Muon::MuonStationIndex::chName
const std::string & chName(ChIndex index)
convert ChIndex into a string
Definition: MuonStationIndex.cxx:119
MuonHough::MuonPhiLayerHough::Maximum
Definition: MuonPhiLayerHough.h:23
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
Muon::MuonLayerHoughTool::Road::maxima
MaximumVec maxima
Definition: MuonLayerHoughTool.h:84
MuonHough::MuonLayerHough::Maximum::triggerConfirmed
int triggerConfirmed
Definition: MuonLayerHough.h:80
Muon::TgcClusterObj3D::phiCluster
HitList phiCluster
Definition: TgcHitClustering.h:27
MuonHough::MuonPhiLayerHough::Maximum::hits
PhiHitVec hits
Definition: MuonPhiLayerHough.h:35
Muon::MuonLayerHoughTool::initializeSectorMapping
void initializeSectorMapping(const EventContext &ctx) const
Definition: MuonLayerHoughTool.cxx:1868
MuonPatternCombination.h
Muon::MuonPrepDataCollection::identify
virtual Identifier identify() const override final
Muon::MuonLayerHoughTool::insertHash
void insertHash(const IdentifierHash &hash, const Identifier &id) const
Definition: MuonLayerHoughTool.cxx:1857
MuonGM::TgcReadoutElement
A TgcReadoutElement corresponds to a single TGC chamber; therefore typically a TGC station contains s...
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/TgcReadoutElement.h:42
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
Muon::MuonStationIndex::ChIndex::BML
@ BML
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
Trk::TrkDetElementBase::surface
virtual const Surface & surface() const =0
Return surface associated with this detector element.
Muon::MuonLayerHoughTool::CollectionsPerSector
Definition: MuonLayerHoughTool.h:53
Muon::MuonLayerHoughTool::createPatternCombinations
void createPatternCombinations(std::vector< MaximumVec > &maxima, MuonPatternCombinationCollection &patternCombis) const
Muon::MuonStationIndex::ChIndex::BOS
@ BOS
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Muon::MuonLayerHoughTool::mergePhiMaxima
void mergePhiMaxima(Road &road) const
Definition: MuonLayerHoughTool.cxx:416
MuonHough::HitDebugInfo::isEtaPhi
bool isEtaPhi
number of layers in the cluster
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:46
MuonHough::MuonDetectorHough::phiHough
MuonPhiLayerHough & phiHough(DetRegIdx region)
access phi transform
Definition: MuonRegionHough.h:71
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Muon::MuonStationIndex::technologyName
const std::string & technologyName(TechnologyIndex index)
convert LayerIndex into a string
Definition: MuonStationIndex.cxx:169
MuonGM::TgcReadoutElement::nStrips
int nStrips(int gasGap) const
Returns the number of strips in a given gas gap.
Muon::MuonPrepDataCollection
Template to hold collections of MuonPrepRawData objects.
Definition: MuonPrepDataCollection.h:46
DataVector< Muon::MuonSegmentCombination >
Muon::MuonLayerHoughTool::Road::mergedPhiMaxima
std::vector< MuonHough::MuonPhiLayerHough::Maximum > mergedPhiMaxima
Definition: MuonLayerHoughTool.h:89
Muon::MuonStationIndex::ChIndex::EES
@ EES
MuonGM::MuonPadDesign::channelWidth
double channelWidth(const Amg::Vector2D &pos, bool measPhi, bool preciseMeas=false) const
calculate local channel width
Definition: MuonPadDesign.h:142
Muon::MuonLayerHoughTool::MaximumVec
HoughDataPerSec::MaximumVec MaximumVec
Definition: MuonLayerHoughTool.h:63
MuonHough::MuonPhiLayerHough::Maximum::max
float max
Definition: MuonPhiLayerHough.h:26
Muon::MuonStationIndex::ChIndex::BEE
@ BEE
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition: MuonStationIndex.h:23
Trk::PrepRawData
Definition: PrepRawData.h:62
MuonHough::MuonPhiLayerHough::Maximum::sector
int sector
Definition: MuonPhiLayerHough.h:33
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
Muon::MuonLayerHoughTool::State::truthHits
std::set< Identifier > truthHits
Definition: MuonLayerHoughTool.h:126
Muon::MuonStationIndex::ChIndex::BMS
@ BMS
Muon::MuonLayerHoughTool::State
Definition: MuonLayerHoughTool.h:123
MuonHough::MuonPhiLayerHough::Maximum::binpos
int binpos
Definition: MuonPhiLayerHough.h:29
MuonHough::PhiHit::prd
const Trk::PrepRawData * prd
access to assiciated hit, either the prd or the tgc pointer is set in athena
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:126
MuonPadDesign.h
Muon::MuonLayerHoughTool::Road
Definition: MuonLayerHoughTool.h:70
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
Muon::MuonLayerHoughTool::associateMaximaInNeighbouringSectors
void associateMaximaInNeighbouringSectors(HoughDataPerSector &houghData, std::vector< HoughDataPerSector > &houghDataPerSectorVec) const
Definition: MuonLayerHoughTool.cxx:780
Muon::MuonStationIndex::layerName
const std::string & layerName(LayerIndex index)
convert LayerIndex into a string
Definition: MuonStationIndex.cxx:153
Muon::MuonLayerHoughTool::Road::neighbouringRegion
DetRegIdx neighbouringRegion
Definition: MuonLayerHoughTool.h:76
Trk::PrepRawData::localPosition
const Amg::Vector2D & localPosition() const
return the local position reference
MuonHough::MuonLayerHough::Maximum::bintheta
int bintheta
Definition: MuonLayerHough.h:79
MuonHough::MuonLayerHough::Maximum::pos
float pos
Definition: MuonLayerHough.h:68
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
MuonHough::MuonDetectorHough
class managing all Hough transforms in the detector
Definition: MuonRegionHough.h:66
checkTriggerxAOD.found
found
Definition: checkTriggerxAOD.py:328
Muon::MdtPrepData
Class to represent measurements from the Monitored Drift Tubes.
Definition: MdtPrepData.h:33
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Muon::MuonLayerHoughTool::associatePhiMaxima
void associatePhiMaxima(Road &road, PhiMaximumVec &phiMaxima) const
Definition: MuonLayerHoughTool.cxx:652
Muon::MuonStationIndex::isBarrel
bool isBarrel(const ChIndex index)
Returns true if the chamber index points to a barrel chamber.
Muon::TgcClusterObj3D
Definition: TgcHitClustering.h:19
Muon::MuonStationIndex::regionName
const std::string & regionName(DetectorRegionIndex index)
convert DetectorRegionIndex into a string
Definition: MuonStationIndex.cxx:138
MuonDetectorManager.h
Muon::MuonLayerHoughTool::printTruthSummary
void printTruthSummary(std::set< Identifier > &truth, std::set< Identifier > &found) const
Definition: MuonLayerHoughTool.cxx:1994
Muon::HoughDataPerSec
Definition: HoughDataPerSec.h:20
Muon::MuonLayerHoughTool::matchTruth
void matchTruth(std::set< Identifier > &truthHits, const PRD_MultiTruthCollection &truthCol, const Identifier &id, MuonHough::HitDebugInfo &debug) const
Definition: MuonLayerHoughTool.cxx:1462
Muon::MuonLayerHoughTool::find
virtual std::pair< std::unique_ptr< MuonPatternCombinationCollection >, std::unique_ptr< HoughDataPerSectorVec > > find(const MdtPrepDataContainer *mdtCont, const CscPrepDataContainer *cscCols, const TgcPrepDataContainer *tgcCont, const RpcPrepDataContainer *rpcCont, const sTgcPrepDataContainer *stgcCont, const MMPrepDataContainer *mmCont, const EventContext &ctx) const override
Definition: MuonLayerHoughTool.cxx:147
Muon::MuonLayerHoughTool::Road::phiMaxima
PhiMaximumVec phiMaxima
Definition: MuonLayerHoughTool.h:85
MuonHough::MuonLayerHough::Maximum::hits
HitVec hits
Definition: MuonLayerHough.h:81
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
CxxUtils::sincos::sn
double sn
Definition: sincos.h:51
MuonGM::MuonChannelDesign::channelLength
double channelLength(int channel) const
STRIPS ONLY: calculate channel length for a given strip number.
Definition: MuonChannelDesign.h:386
MuonChannelDesign.h
Muon::MuonLayerHoughTool::fill
void fill(const EventContext &ctx, std::set< Identifier > &truthHits, const MdtPrepDataCollection &mdts, HitVec &hits) const
Definition: MuonLayerHoughTool.cxx:1481
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
IdentifiableContainerMT::indexFindPtr
virtual const T * indexFindPtr(IdentifierHash hashId) const override final
return pointer on the found entry or null if out of range using hashed index - fast version,...
Definition: IdentifiableContainerMT.h:289
Muon::MuonLayerHoughTool::fillHitsPerSector
void fillHitsPerSector(const EventContext &ctx, State &state, const int sector, const CollectionsPerSector &hashes, const MdtPrepDataContainer *mdtCont, const CscPrepDataContainer *cscCont, const TgcPrepDataContainer *tgcCont, const RpcPrepDataContainer *rpcCont, const sTgcPrepDataContainer *stgcCont, const MMPrepDataContainer *mmCont) const
Definition: MuonLayerHoughTool.cxx:1421
Muon::MuonLayerHoughTool::PhiHitVec
HoughDataPerSec::PhiHitVec PhiHitVec
Definition: MuonLayerHoughTool.h:61
MuonGM::MuonChannelDesign
Definition: MuonChannelDesign.h:24
Muon::TgcHitClusteringObj
Definition: TgcHitClustering.h:47
Muon::MuonStationIndex::LayerIndex::LayerIndexMax
@ LayerIndexMax
BEE.
MuonHough::Hit::tgc
const Muon::TgcClusterObj3D * tgc
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:89
python.selector.AtlRunQuerySelectorLhcOlc.selector
selector
Definition: AtlRunQuerySelectorLhcOlc.py:610
Muon::MuonStationIndex::DetectorRegionIndex
DetectorRegionIndex
enum to classify the different layers in the muon spectrometer
Definition: MuonStationIndex.h:47
y
#define y
Muon::HoughDataPerSec::maxAssociationMap
MaximumAssociationMap maxAssociationMap
Definition: HoughDataPerSec.h:60
CaloCondBlobAlgs_fillNoiseFromASCII.hash
dictionary hash
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:108
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
Definition: GeoPrimitivesHelpers.h:347
lumiFormat.fill
fill
Definition: lumiFormat.py:104
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Muon::MuonLayerHoughTool::CollectionsPerSector::technologyRegionHashVecs
TechnologyRegionHashVec technologyRegionHashVecs
Definition: MuonLayerHoughTool.h:55
DeMoScan.first
bool first
Definition: DeMoScan.py:534
DEBUG
#define DEBUG
Definition: page_access.h:11
Muon::MuonStationIndex::TechnologyIndex::STGC
@ STGC
Muon::TgcPrepData
Class to represent TGC measurements.
Definition: TgcPrepData.h:32
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:39
Muon::MuonStationIndex::ChIndex::EOL
@ EOL
Muon::MuonLayerHoughTool::buildRoads
void buildRoads(MaximumVec &seedMaxima, MuonHough::MuonDetectorHough &detectorHoughTransforms, std::unique_ptr< HoughDataPerSectorVec > &houghDataPerSectorVec, std::vector< Road > &roads) const
Definition: MuonLayerHoughTool.cxx:300
Muon::HoughDataPerSec::sector
int sector
Definition: HoughDataPerSec.h:51
MuonHough::Hit::prd
const Trk::PrepRawData * prd
access to assiciated hit, either the prd or the tgc pointer is set in athena
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:88
Muon::MuonLayerHoughTool::CollectionsPerSector::sector
int sector
Definition: MuonLayerHoughTool.h:54
MuonHough::extrapolate
float extrapolate(const MuonLayerHough::Maximum &ref, const MuonLayerHough::Maximum &ex, bool doparabolic=false)
Definition: MuonLayerHough.cxx:521
Muon::MuonStationIndex::toLayerIndex
LayerIndex toLayerIndex(ChIndex index)
convert ChIndex into LayerIndex
Muon::MuonStationIndex::ChIndex
ChIndex
enum to classify the different chamber layers in the muon spectrometer
Definition: MuonStationIndex.h:15
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:13
MuonHough::MuonPhiLayerHough::setDebug
void setDebug(bool d)
Definition: MuonPhiLayerHough.h:44
Muon::MuonStationIndex::ChIndex::ChUnknown
@ ChUnknown
Muon::MuonStationIndex::ChIndex::EMS
@ EMS
IdentifierHash
This is a "hash" representation of an Identifier. This encodes a 32 bit index which can be used to lo...
Definition: IdentifierHash.h:25
Muon::MuonPatternCombination
The MuonPatternCombination class provides the means to store the output of the initial global pattern...
Definition: MuonPatternCombination.h:29
TruthParticle.h
MuonHough::MuonLayerHough::Maximum::binposmax
int binposmax
Definition: MuonLayerHough.h:77
set_intersection
Set * set_intersection(Set *set1, Set *set2)
Perform an intersection of two sets.
MuonGM::MuonChannelDesign::channelWidth
double channelWidth() const
calculate local channel width
Definition: MuonChannelDesign.h:394
MuonHough::SortHitsPerLayer
struct to sort the hits
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:180
python.SystemOfUnits.m2
float m2
Definition: SystemOfUnits.py:107
Muon::sTgcPrepData
Class to represent sTgc measurements.
Definition: sTgcPrepData.h:20
Muon::MuonStationIndex::decomposeSectorLayerHash
std::pair< DetectorRegionIndex, LayerIndex > decomposeSectorLayerHash(unsigned int hash)
decompose the hash into Region and Layer
dq_make_web_display.cl
cl
print [x.__class__ for x in toList(dqregion.getSubRegions()) ]
Definition: dq_make_web_display.py:25
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
Muon::MuonLayerHoughTool::Road::neighbouringSector
int neighbouringSector
Definition: MuonLayerHoughTool.h:77
Muon::HoughDataPerSec::hitVec
RegionHitVec hitVec
Definition: HoughDataPerSec.h:52
MuonHough::PhiHit
struct containing all hit information needed for the Hough transform
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:99
Muon::TgcHitClusteringObj::clusters3D
std::vector< TgcClusterObj3D > clusters3D
Definition: TgcHitClustering.h:69
Muon::MuonStationIndex::StIndex::EI
@ EI
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
Trk::Surface::localToGlobal
virtual void localToGlobal(const Amg::Vector2D &locp, const Amg::Vector3D &mom, Amg::Vector3D &glob) const =0
Specified by each surface type: LocalToGlobal method without dynamic memory allocation.
MuonHough::MuonPhiLayerHough::findMaximum
bool findMaximum(Maximum &maximum, float maxval) const
Definition: MuonPhiLayerHough.cxx:190
Muon::TgcClusterObj3D::etaCluster
HitList etaCluster
Definition: TgcHitClustering.h:26
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
Muon::MuonLayerHoughTool::State::houghDataPerSectorVec
std::unique_ptr< HoughDataPerSectorVec > houghDataPerSectorVec
Definition: MuonLayerHoughTool.h:125
ymax
double ymax
Definition: listroot.cxx:64
Muon::MuonStationIndex::toInt
constexpr int toInt(const EnumType enumVal)
Definition: MuonStationIndex.h:61
Trk::PrepRawData::detectorElement
virtual const TrkDetElementBase * detectorElement() const =0
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Muon::HoughDataPerSec::maxVec
RegionMaximumVec maxVec
Definition: HoughDataPerSec.h:54
Muon::MuonStationIndex::ChIndex::BOL
@ BOL
Identifier
Definition: IdentifierFieldParser.cxx:14
Muon::MuonLayerHoughTool::State::outputTruthHits
std::set< Identifier > outputTruthHits
Definition: MuonLayerHoughTool.h:128