ATLAS Offline Software
Loading...
Searching...
No Matches
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
5#include "MuonHoughPatternTools/MuonLayerHoughTool.h"
6
7#include "AtlasHepMC/GenEvent.h"
8#include "CxxUtils/sincos.h"
9#include "GaudiKernel/ConcurrencyFlags.h"
10#include "MuonPattern/MuonPatternChamberIntersect.h"
11#include "MuonPattern/MuonPatternCombination.h"
12#include "MuonReadoutGeometry/MuonChannelDesign.h"
13#include "MuonReadoutGeometry/MuonDetectorManager.h"
14#include "StoreGate/ReadCondHandle.h"
15#include "MuonReadoutGeometry/MuonPadDesign.h"
16#include "MuonReadoutGeometry/sTgcReadoutElement.h"
17#include "xAODMuon/MuonSegmentContainer.h"
18#include "xAODTruth/TruthParticle.h"
19#include "xAODTruth/TruthParticleContainer.h"
20namespace Muon {
21 using namespace MuonStationIndex;
22
23 StatusCode MuonLayerHoughTool::initialize() {
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
191 MuonHough::MuonLayerHough& hough =
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;
227 LayerIndex layer = toLayerIndex(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;
318 MuonStationIndex::LayerIndex layer = MuonStationIndex::toLayerIndex(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;
383 MuonStationIndex::LayerIndex layer = MuonStationIndex::toLayerIndex(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
416 void MuonLayerHoughTool::mergePhiMaxima(MuonLayerHoughTool::Road& road) const {
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
500 void MuonLayerHoughTool::extendSeed(MuonHough::MuonDetectorHough& detectorHoughTransforms,
501 MuonLayerHoughTool::Road& road,
502 MuonLayerHoughTool::HoughDataPerSector& sectorData) const { // const {
503 if (!road.seed) return;
504
505 RegionMaximumVec& maxVec = sectorData.maxVec;
506
507 // gather locality information on seed
508 MuonHough::MuonLayerHough::Maximum& seed = *road.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) {
514 MuonStationIndex::LayerIndex layer = static_cast<MuonStationIndex::LayerIndex>(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) {
566 MuonStationIndex::LayerIndex layer = static_cast<MuonStationIndex::LayerIndex>(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
652 void MuonLayerHoughTool::associatePhiMaxima(MuonLayerHoughTool::Road& road, MuonLayerHoughTool::PhiMaximumVec& phiMaxima) const {
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
780 void MuonLayerHoughTool::associateMaximaInNeighbouringSectors(
781 MuonLayerHoughTool::HoughDataPerSector& houghData,
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
855 void MuonLayerHoughTool::associateMaximaToPhiMaxima(
856 DetRegIdx region, MuonLayerHoughTool::HoughDataPerSector& houghData,
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) {
911 MuonStationIndex::LayerIndex layer = static_cast<MuonStationIndex::LayerIndex>(lay);
912 unsigned int layerHash = MuonStationIndex::sectorLayerHash(region, layer);
913 MaximumVec& maxima = houghData.maxVec[layerHash];
914 if (maxima.empty()) continue;
915 MuonStationIndex::StIndex stIndex = MuonStationIndex::toStationIndex(region, layer);
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) {
1027 MuonStationIndex::LayerIndex layer = static_cast<MuonStationIndex::LayerIndex>(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
1046 void MuonLayerHoughTool::createPatternCombinations(MuonLayerHoughTool::RegionMaximumVec& maxima,
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
1108 void MuonLayerHoughTool::createPatternCombinations(
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);
1210 ATH_MSG_DEBUG(MuonStationIndex::chName(chIndex)
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);
1243 Amg::Vector3D gpos{Amg::Vector3D::Zero()};
1244 Amg::Vector3D gdir{Amg::Vector3D::Zero()};
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
1277 bool MuonLayerHoughTool::findMaxima(MaximumVec& seedMaxima,
1278 MuonHough::MuonLayerHough& hough, MuonLayerHoughTool::HitVec& hits,
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;
1296 MuonHough::MuonLayerHoughSelector selector;
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) {
1309 MuonHough::MuonLayerHough::Maximum maximum;
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
1351 bool MuonLayerHoughTool::findMaxima(MuonHough::MuonPhiLayerHough& hough, MuonLayerHoughTool::PhiHitVec& hits,
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) {
1361 MuonHough::MuonPhiLayerHough::Maximum maximum;
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,
1482 MuonLayerHoughTool::HitVec& hits) const {
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 << " "
1516 << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
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 << " "
1538 << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
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,
1570 MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits) const {
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 << " "
1586 << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
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,
1620 MuonLayerHoughTool::HitVec& hits) const {
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 << " "
1675 << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
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
1839 MuonHough::HitDebugInfo* phiDebug = new MuonHough::HitDebugInfo(*debug);
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 << " "
1852 << MuonStationIndex::regionName(region) << " " << MuonStationIndex::layerName(layer)
1853 << " -> etaHits: " << hits.size() << " phiHits: " << phiHits.size()
1854 << " sectors: " << sectors.size());
1855 }
1856
1857 void MuonLayerHoughTool::insertHash(const IdentifierHash& hash, const Identifier& id) const{
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 "
1887 << MuonStationIndex::sectorLayerHashMax());
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) {
1894 IdentifierHash hash;
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) {
1914 IdentifierHash hash;
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) {
1924 IdentifierHash hash;
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);
1941 IdentifierHash hash;
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 =
1970 MuonStationIndex::decomposeSectorLayerHash(hash);
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)
1984 << MuonStationIndex::technologyName(static_cast<MuonStationIndex::TechnologyIndex>(tech))
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());
2000 std::vector<Identifier>::iterator pos =
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
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
vec
std::vector< size_t > vec
Definition CombinationsGeneratorTest.cxx:9
debug
const bool debug
Definition MakeUncertaintyPlots.cxx:53
MuonPatternCombinationCollection
DataVector< Muon::MuonPatternCombination > MuonPatternCombinationCollection
This typedef represents a collection of MuonPatternCombination objects.
Definition MuonPatternCombinationCollection.h:17
MuonSegmentCombinationCollection
DataVector< Muon::MuonSegmentCombination > MuonSegmentCombinationCollection
This typedef represents a collection of MuonSegmentCombination objects.
Definition MuonSegmentCombinationCollection.h:16
sign
int sign(int a)
Definition TRT_StrawNeighbourSvc.h:108
y
#define y
x
#define x
z
#define z
empty
static const Attributes_t empty
Definition XmlStreamer.cxx:16
max
#define max(a, b)
Definition cfImp.cxx:41
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DataVector< PrepDataT >::end
const_iterator end() const noexcept
DataVector< PrepDataT >::begin
const_iterator begin() const noexcept
DataVector< PrepDataT >::size
size_type size() const noexcept
DataVector< PrepDataT >::empty
bool empty() const noexcept
IdentifiableContainerMT::empty
bool empty() const
return true if container is empty
Definition IdentifiableContainerMT.h:245
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:290
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
MuonGM::TgcReadoutElement
A TgcReadoutElement corresponds to a single TGC chamber; therefore typically a TGC station contains s...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/TgcReadoutElement.h:42
MuonGM::TgcReadoutElement::channelPos
Amg::Vector3D channelPos(const Identifier &id) const
Returns the position of the active channel (wireGang or strip)
MuonGM::TgcReadoutElement::nStrips
int nStrips(int gasGap) const
Returns the number of strips in a given gas gap.
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
MuonHough::Hit
struct containing all hit information needed for the Hough transform
Definition MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:60
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
MuonHough::MuonDetectorHough
class managing all Hough transforms in the detector
Definition MuonRegionHough.h:66
MuonHough::MuonDetectorHough::phiHough
MuonPhiLayerHough & phiHough(DetRegIdx region)
access phi transform
Definition MuonRegionHough.h:71
Muon::CscPrepData
Class representing clusters from the CSC.
Definition CscPrepData.h:39
Muon::MMPrepData
Class to represent MM measurements.
Definition MMPrepData.h:22
Muon::MdtPrepData
Class to represent measurements from the Monitored Drift Tubes.
Definition MdtPrepData.h:33
Muon::MuonLayerHoughTool::Road::add
void add(std::shared_ptr< MuonHough::MuonLayerHough::Maximum > max)
Definition MuonLayerHoughTool.h:79
Muon::MuonLayerHoughTool::Road::mergedPhiMaxima
std::vector< MuonHough::MuonPhiLayerHough::Maximum > mergedPhiMaxima
Definition MuonLayerHoughTool.h:89
Muon::MuonLayerHoughTool::Road::maximumSet
std::set< std::shared_ptr< MuonHough::MuonLayerHough::Maximum > > maximumSet
Definition MuonLayerHoughTool.h:87
Muon::MuonLayerHoughTool::Road::seed
std::shared_ptr< MuonHough::MuonLayerHough::Maximum > seed
Definition MuonLayerHoughTool.h:78
Muon::MuonLayerHoughTool::m_useSeeds
Gaudi::Property< bool > m_useSeeds
Definition MuonLayerHoughTool.h:187
Muon::MuonLayerHoughTool::printTruthSummary
void printTruthSummary(std::set< Identifier > &truth, std::set< Identifier > &found) const
Definition MuonLayerHoughTool.cxx:1994
Muon::MuonLayerHoughTool::RegionMaximumVec
HoughDataPerSec::RegionMaximumVec RegionMaximumVec
Definition MuonLayerHoughTool.h:66
Muon::MuonLayerHoughTool::m_selectorsLoose
std::vector< MuonHough::MuonLayerHoughSelector > m_selectorsLoose
Definition MuonLayerHoughTool.h:196
Muon::MuonLayerHoughTool::m_selectors
std::vector< MuonHough::MuonLayerHoughSelector > m_selectors
Definition MuonLayerHoughTool.h:195
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::MuonLayerHoughTool::extendSeed
void extendSeed(MuonHough::MuonDetectorHough &detectorHoughTransforms, Road &road, HoughDataPerSector &sectorData) const
Definition MuonLayerHoughTool.cxx:500
Muon::MuonLayerHoughTool::mergePhiMaxima
void mergePhiMaxima(Road &road) const
Definition MuonLayerHoughTool.cxx:416
Muon::MuonLayerHoughTool::HitVec
HoughDataPerSec::HitVec HitVec
Definition MuonLayerHoughTool.h:59
Muon::MuonLayerHoughTool::m_truthNames
SG::ReadHandleKeyArray< PRD_MultiTruthCollection > m_truthNames
Definition MuonLayerHoughTool.h:198
Muon::MuonLayerHoughTool::findMaxima
bool findMaxima(MaximumVec &seedMaxima, MuonHough::MuonLayerHough &hough, HitVec &hits, MaximumVec &maxima) const
Definition MuonLayerHoughTool.cxx:1277
Muon::MuonLayerHoughTool::HoughDataPerSector
HoughDataPerSec HoughDataPerSector
Definition MuonLayerHoughTool.h:68
Muon::MuonLayerHoughTool::getSectors
void getSectors(const Amg::Vector3D &pos, std::vector< int > &sectors) const
Definition MuonLayerHoughTool.h:228
Muon::MuonLayerHoughTool::m_idHelperSvc
ServiceHandle< IMuonIdHelperSvc > m_idHelperSvc
Definition MuonLayerHoughTool.h:192
Muon::MuonLayerHoughTool::m_sectorMapping
MuonSectorMapping m_sectorMapping
Definition MuonLayerHoughTool.h:211
Muon::MuonLayerHoughTool::m_printer
PublicToolHandle< MuonEDMPrinterTool > m_printer
Definition MuonLayerHoughTool.h:193
Muon::MuonLayerHoughTool::analyse
std::pair< std::unique_ptr< MuonPatternCombinationCollection >, std::unique_ptr< HoughDataPerSectorVec > > analyse(State &state) const
Definition MuonLayerHoughTool.cxx:170
Muon::MuonLayerHoughTool::m_requireTriggerConfirmationNSW
Gaudi::Property< bool > m_requireTriggerConfirmationNSW
Definition MuonLayerHoughTool.h:201
Muon::MuonLayerHoughTool::sublay
int sublay(const Identifier &id, float z=0) const
Definition MuonLayerHoughTool.h:245
Muon::MuonLayerHoughTool::fill
void fill(const EventContext &ctx, std::set< Identifier > &truthHits, const MdtPrepDataCollection &mdts, HitVec &hits) const
Definition MuonLayerHoughTool.cxx:1481
Muon::MuonLayerHoughTool::m_techToTruthNameIdx
std::map< unsigned int, unsigned int > m_techToTruthNameIdx
Definition MuonLayerHoughTool.h:209
Muon::MuonLayerHoughTool::DetRegIdx
MuonStationIndex::DetectorRegionIndex DetRegIdx
Definition MuonLayerHoughTool.h:51
Muon::MuonLayerHoughTool::createPatternCombinations
void createPatternCombinations(std::vector< MaximumVec > &maxima, MuonPatternCombinationCollection &patternCombis) const
Muon::MuonLayerHoughTool::rCor
double rCor(const Amg::Vector3D &pos, const Identifier &id) const
Definition MuonLayerHoughTool.h:232
Muon::MuonLayerHoughTool::m_muonManagerKey
SG::ReadCondHandleKey< MuonGM::MuonDetectorManager > m_muonManagerKey
Definition MuonLayerHoughTool.h:189
Muon::MuonLayerHoughTool::PhiMaximumVec
HoughDataPerSec::PhiMaximumVec PhiMaximumVec
Definition MuonLayerHoughTool.h:64
Muon::MuonLayerHoughTool::TechnologyRegionHashVec
std::vector< RegionHashVec > TechnologyRegionHashVec
Definition MuonLayerHoughTool.h:50
Muon::MuonLayerHoughTool::insertHash
void insertHash(const IdentifierHash &hash, const Identifier &id) const
Definition MuonLayerHoughTool.cxx:1857
Muon::MuonLayerHoughTool::PhiHitVec
HoughDataPerSec::PhiHitVec PhiHitVec
Definition MuonLayerHoughTool.h:61
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::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::m_extrapolationDistance
Gaudi::Property< float > m_extrapolationDistance
Definition MuonLayerHoughTool.h:206
Muon::MuonLayerHoughTool::associatePhiMaxima
void associatePhiMaxima(Road &road, PhiMaximumVec &phiMaxima) const
Definition MuonLayerHoughTool.cxx:652
Muon::MuonLayerHoughTool::initialize
virtual StatusCode initialize() override
Definition MuonLayerHoughTool.cxx:23
Muon::MuonLayerHoughTool::m_addSectors
Gaudi::Property< bool > m_addSectors
Definition MuonLayerHoughTool.h:207
Muon::MuonLayerHoughTool::m_debugHough
Gaudi::Property< bool > m_debugHough
Definition MuonLayerHoughTool.h:203
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::associateMaximaInNeighbouringSectors
void associateMaximaInNeighbouringSectors(HoughDataPerSector &houghData, std::vector< HoughDataPerSector > &houghDataPerSectorVec) const
Definition MuonLayerHoughTool.cxx:780
Muon::MuonLayerHoughTool::m_doParabolicExtrapolation
Gaudi::Property< bool > m_doParabolicExtrapolation
Definition MuonLayerHoughTool.h:204
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::MuonLayerHoughTool::HashVec
std::vector< IdentifierHash > HashVec
Definition MuonLayerHoughTool.h:48
Muon::MuonLayerHoughTool::initializeSectorMapping
void initializeSectorMapping(const EventContext &ctx) const
Definition MuonLayerHoughTool.cxx:1868
Muon::MuonLayerHoughTool::MaximumVec
HoughDataPerSec::MaximumVec MaximumVec
Definition MuonLayerHoughTool.h:63
Muon::MuonPatternChamberIntersect
This class holds information needed for the Moore and MoMu pattern recognition for a muon chamber.
Definition MuonPatternChamberIntersect.h:38
Muon::MuonPatternCombination
The MuonPatternCombination class provides the means to store the output of the initial global pattern...
Definition MuonPatternCombination.h:29
Muon::MuonPrepDataCollection::identify
virtual Identifier identify() const override final
Definition MuonPrepDataCollection.h:39
Muon::RpcPrepData
Class to represent RPC measurements.
Definition RpcPrepData.h:35
Muon::TgcPrepData
Class to represent TGC measurements.
Definition TgcPrepData.h:32
Muon::sTgcPrepData
Class to represent sTgc measurements.
Definition sTgcPrepData.h:20
PRD_MultiTruthCollection
A PRD is mapped onto all contributing particles.
Definition PRD_MultiTruthCollection.h:24
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...
Trk::PrepRawData::localPosition
const Amg::Vector2D & localPosition() const
return the local position reference
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
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.
Trk::TrkDetElementBase::surface
virtual const Surface & surface() const =0
Return surface associated with this detector element.
intersection
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Definition compareFlatTrees.cxx:25
r
int r
Definition globals.cxx:22
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146
ymin
double ymin
Definition listroot.cxx:63
ymax
double ymax
Definition listroot.cxx:64
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
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition GeoPrimitives.h:48
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
HepMC::uniqueID
int uniqueID(const T &p)
Definition MagicNumbers.h:117
MuonHough::PhiHitVec
std::vector< std::shared_ptr< MuonHough::PhiHit > > PhiHitVec
Definition MuonPhiLayerHough.h:20
MuonHough::extrapolate
float extrapolate(const MuonLayerHough::Maximum &ref, const MuonLayerHough::Maximum &ex, bool doparabolic=false)
Definition MuonLayerHough.cxx:521
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition MuonStationIndex.h:23
Muon::MuonStationIndex::layerName
const std::string & layerName(LayerIndex index)
convert LayerIndex into a string
Definition MuonStationIndex.cxx:153
Muon::MuonStationIndex::chIndex
ChIndex chIndex(const std::string &index)
convert ChIndex name string to enum
Definition MuonStationIndex.cxx:11
Muon::MuonStationIndex::TechnologyIndex
TechnologyIndex
enum to classify the different layers in the muon spectrometer
Definition MuonStationIndex.h:54
Muon::MuonStationIndex::DetectorRegionIndex
DetectorRegionIndex
enum to classify the different layers in the muon spectrometer
Definition MuonStationIndex.h:47
Muon::MuonStationIndex::toStationIndex
StIndex toStationIndex(ChIndex index)
convert ChIndex into StIndex
Muon::MuonStationIndex::numberOfSectors
constexpr unsigned numberOfSectors()
return total number of sectors
Definition MuonStationIndex.h:118
Muon::MuonStationIndex::decomposeSectorLayerHash
std::pair< DetectorRegionIndex, LayerIndex > decomposeSectorLayerHash(unsigned int hash)
decompose the hash into Region and Layer
Muon::MuonStationIndex::toInt
constexpr int toInt(const EnumType enumVal)
Definition MuonStationIndex.h:61
Muon::MuonStationIndex::sectorLayerHashMax
constexpr unsigned int sectorLayerHashMax()
maximum create a hash out of region and layer
Definition MuonStationIndex.h:110
Muon::MuonStationIndex::isBarrel
bool isBarrel(const ChIndex index)
Returns true if the chamber index points to a barrel chamber.
Muon::MuonStationIndex::sectorLayerHash
unsigned int sectorLayerHash(DetectorRegionIndex detectorRegionIndex, LayerIndex layerIndex)
create a hash out of region and layer
Muon::MuonStationIndex::stName
const std::string & stName(StIndex index)
convert StIndex into a string
Definition MuonStationIndex.cxx:104
Muon::MuonStationIndex::technologyName
const std::string & technologyName(TechnologyIndex index)
convert LayerIndex into a string
Definition MuonStationIndex.cxx:169
Muon::MuonStationIndex::chName
const std::string & chName(ChIndex index)
convert ChIndex into a string
Definition MuonStationIndex.cxx:119
Muon::MuonStationIndex::regionName
const std::string & regionName(DetectorRegionIndex index)
convert DetectorRegionIndex into a string
Definition MuonStationIndex.cxx:138
Muon::MuonStationIndex::toLayerIndex
LayerIndex toLayerIndex(ChIndex index)
convert ChIndex into LayerIndex
Muon::MuonStationIndex::LayerIndex
LayerIndex
enum to classify the different layers in the muon spectrometer
Definition MuonStationIndex.h:38
Muon::MuonStationIndex::ChIndex
ChIndex
enum to classify the different chamber layers in the muon spectrometer
Definition MuonStationIndex.h:15
Muon
NRpcCablingAlg reads raw condition data and writes derived condition data to the condition store.
Definition TrackSystemController.h:45
Muon::CscPrepDataCollection
MuonPrepDataCollection< CscPrepData > CscPrepDataCollection
Definition MuonPrepDataCollection.h:109
Muon::TgcPrepDataCollection
MuonPrepDataCollection< TgcPrepData > TgcPrepDataCollection
Definition MuonPrepDataCollection.h:108
Muon::RpcPrepDataContainer
MuonPrepDataContainerT< RpcPrepData > RpcPrepDataContainer
Definition MuonPrepDataContainer.h:96
Muon::MMPrepDataCollection
MuonPrepDataCollection< MMPrepData > MMPrepDataCollection
Definition MuonPrepDataCollection.h:113
Muon::TgcPrepDataContainer
MuonPrepDataContainerT< TgcPrepData > TgcPrepDataContainer
Definition MuonPrepDataContainer.h:97
Muon::MdtPrepDataContainer
MuonPrepDataContainerT< MdtPrepData > MdtPrepDataContainer
Definition MuonPrepDataContainer.h:95
Muon::sTgcPrepDataContainer
MuonPrepDataContainerT< sTgcPrepData > sTgcPrepDataContainer
Definition MuonPrepDataContainer.h:101
Muon::MdtStatusDriftTime
@ MdtStatusDriftTime
The tube produced a vaild measurement.
Definition MdtDriftCircleStatus.h:34
Muon::MMPrepDataContainer
MuonPrepDataContainerT< MMPrepData > MMPrepDataContainer
Definition MuonPrepDataContainer.h:100
Muon::MdtPrepDataCollection
MuonPrepDataCollection< MdtPrepData > MdtPrepDataCollection
Definition MuonPrepDataCollection.h:106
Muon::CscPrepDataContainer
MuonPrepDataContainerT< CscPrepData > CscPrepDataContainer
Definition MuonPrepDataContainer.h:98
Muon::RpcPrepDataCollection
MuonPrepDataCollection< RpcPrepData > RpcPrepDataCollection
Definition MuonPrepDataCollection.h:107
Muon::sTgcPrepDataCollection
MuonPrepDataCollection< sTgcPrepData > sTgcPrepDataCollection
Definition MuonPrepDataCollection.h:112
Trk::locR
@ locR
Definition ParamDefs.h:44
index
Definition index.py:1
std::stable_sort
void stable_sort(DataModel_detail::iterator< DVL > beg, DataModel_detail::iterator< DVL > end)
Specialization of stable_sort for DataVector/List.
Definition DVL_algorithms.h:644
xAOD::P4Helpers::deltaPhi
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
Definition xAODP4Helpers.h:69
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition sincos.h:39
MuonGM::MuonChannelDesign::channelWidth
double channelWidth() const
calculate local channel width
Definition MuonChannelDesign.h:391
MuonGM::MuonChannelDesign::channelLength
double channelLength(int channel) const
STRIPS ONLY: calculate channel length for a given strip number.
Definition MuonChannelDesign.h:383
MuonGM::MuonPadDesign
Parameters defining the design of the readout sTGC pads.
Definition MuonPadDesign.h:37
MuonGM::MuonPadDesign::channelWidth
double channelWidth(const Amg::Vector2D &pos, bool measPhi, bool preciseMeas=false) const
calculate local channel width
Definition MuonPadDesign.h:139
MuonHough::MuonLayerHough::Maximum
struct representing the maximum in the hough space
Definition MuonLayerHough.h:64
MuonHough::MuonLayerHough::m_descriptor
RegionDescriptor m_descriptor
Definition MuonLayerHough.h:187
MuonHough::MuonPhiLayerHough::fillLayer2
void fillLayer2(const PhiHitVec &hits, bool substract=false) const
Definition MuonPhiLayerHough.cxx:29
MuonHough::MuonPhiLayerHough::findMaximum
bool findMaximum(Maximum &maximum, float maxval) const
Definition MuonPhiLayerHough.cxx:190
MuonHough::MuonPhiLayerHough::associateHitsToMaximum
void associateHitsToMaximum(Maximum &maximum, const PhiHitVec &hits) const
Definition MuonPhiLayerHough.cxx:246
MuonHough::PhiHit
struct containing all hit information needed for the Hough transform
Definition MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:99
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
Muon::HoughDataPerSec::nlayersWithMaxima
std::array< int, detRegMax > nlayersWithMaxima
Definition HoughDataPerSec.h:56
Muon::HoughDataPerSec::nphimaxHitsInRegion
std::array< int, detRegMax > nphimaxHitsInRegion
Definition HoughDataPerSec.h:59
Muon::HoughDataPerSec::phiMaxVec
RegionPhiMaximumVec phiMaxVec
Definition HoughDataPerSec.h:55
Muon::HoughDataPerSec::maxAssociationMap
MaximumAssociationMap maxAssociationMap
Definition HoughDataPerSec.h:60
Muon::HoughDataPerSec::phiHitVec
RegionPhiHitVec phiHitVec
Definition HoughDataPerSec.h:53
Muon::HoughDataPerSec::nmaxHitsInRegion
std::array< int, detRegMax > nmaxHitsInRegion
Definition HoughDataPerSec.h:58
Muon::HoughDataPerSec::associatedToOtherSector
std::set< MuonHough::MuonLayerHough::Maximum * > associatedToOtherSector
Definition HoughDataPerSec.h:62
Muon::HoughDataPerSec::maxVec
RegionMaximumVec maxVec
Definition HoughDataPerSec.h:54
Muon::HoughDataPerSec::nphilayersWithMaxima
std::array< int, detRegMax > nphilayersWithMaxima
Definition HoughDataPerSec.h:57
Muon::MuonLayerHoughTool::State::houghDataPerSectorVec
std::unique_ptr< HoughDataPerSectorVec > houghDataPerSectorVec
Definition MuonLayerHoughTool.h:125
Muon::MuonLayerHoughTool::State::foundTruthHits
std::set< Identifier > foundTruthHits
Definition MuonLayerHoughTool.h:127
Muon::MuonLayerHoughTool::State::outputTruthHits
std::set< Identifier > outputTruthHits
Definition MuonLayerHoughTool.h:128
Muon::TgcHitClusteringObj::clusters3D
std::vector< TgcClusterObj3D > clusters3D
Definition TgcHitClustering.h:69
Muon::TgcHitClusteringObj::bestEtaCluster
const HitList & bestEtaCluster() const
Definition TgcHitClustering.cxx:12
Muon::TgcHitClusteringObj::cluster
bool cluster(const HitList &col)
result
Generated on for ATLAS Offline Software by doxygen 1.14.0