ATLAS Offline Software
Loading...
Searching...
No Matches
MuonCombinePatternTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "MuonCombinePatternTools/MuonCombinePatternTool.h"
6
7#include "CxxUtils/sincos.h"
8#include "MuonHoughPatternEvent/MuonHoughPattern.h"
9#include "MuonPattern/MuonPatternChamberIntersect.h"
10#include "MuonPrepRawData/MuonPrepDataContainer.h"
11#include "MuonReadoutGeometry/MuonDetectorManager.h"
12#include "TrkParameters/TrackParameters.h"
13#include "TrkSurfaces/Surface.h"
14#include <cmath> //std::abs, M_PI etc
15#include "GaudiKernel/SystemOfUnits.h"
16#include "EventPrimitives/EventPrimitivesToStringConverter.h"
17namespace{
18 double rotatePhi(double phi, double rotationFraction) {
19 // check whether we rotate to a large value than pi, if so add additional
20 // rotation by -2pi
21 if (phi + rotationFraction * M_PI > M_PI) return phi + (rotationFraction - 2.) * M_PI;
22 return phi + rotationFraction * M_PI;
23 }
24
25 struct Unowned{
26 void operator()(const Muon::MuonPrdPattern*) {}
27 };
28
29 const Amg::Vector3D& globalPos(const Trk::PrepRawData* prd ){
30 if (prd->type(Trk::PrepRawDataType::MdtPrepData))
31 return static_cast<const Muon::MdtPrepData*>(prd)->globalPosition();
32 return static_cast<const Muon::MuonCluster*>(prd)->globalPosition();
33
34 }
35}
36using PrdPatternPair = MuonCombinePatternTool::PrdPatternPair;
37
38MuonCombinePatternTool::MuonCombinePatternTool(const std::string& type, const std::string& name, const IInterface* parent) :
39 AthAlgTool(type, name, parent),
40 m_maximum_xydistance(3500),
41 m_maximum_rzdistance(1500),
42 m_use_cosmics(false),
43 m_splitpatterns(true),
44 m_nodiscarding(true),
45 m_bestphimatch(false),
46 m_flipdirectionforcosmics(false) {
47 declareInterface<IMuonCombinePatternTool>(this);
48 declareProperty("UseCosmics", m_use_cosmics);
49 declareProperty("SplitPatterns", m_splitpatterns);
50 declareProperty("NoDiscarding", m_nodiscarding);
51 declareProperty("BestPhiMatch", m_bestphimatch);
52 declareProperty("FlipDirectionForCosmics", m_flipdirectionforcosmics);
53 declareProperty("UseTightAssociation", m_useTightAssociation = false);
54 declareProperty("MaxSizePhiPatternLooseCuts", m_maxSizePhiPatternLoose = 40);
55 declareProperty("MaxSizeEtaPatternLooseCuts", m_maxSizeEtaPatternLoose = 200);
56}
57
58StatusCode MuonCombinePatternTool::initialize() {
59 ATH_MSG_DEBUG("MuonCombinePatternTool::initialize");
60 ATH_CHECK(m_idHelperSvc.retrieve());
61 ATH_CHECK(m_printer.retrieve());
62 if (!m_use_cosmics) { m_splitpatterns = false; }
63 if (m_use_cosmics) { m_bestphimatch = true; }
64 ATH_MSG_DEBUG(" UseCosmics: " << m_use_cosmics << " Split Patterns: " << m_splitpatterns << " NoDiscarding: " << m_nodiscarding
65 << " BestPhiMatch: " << m_bestphimatch);
66 return StatusCode::SUCCESS;
67}
68
69std::unique_ptr<MuonPrdPatternCollection> MuonCombinePatternTool::combineEtaPhiPatterns(
70 const MuonPrdPatternCollection& phiPatternCollection, const MuonPrdPatternCollection& etaPatternCollection,
72 const EtaPhiHitAssocMap& phiEtaHitAssMap) const {
73 bool myDebug = false;
74
78 std::vector<CandidatePatPair> candidates{};
79
80
81 // strategy
82 // associate eta pattern to phi patterns
83 // are phi hits on a eta and eta hits on phi pattern?
84
85 // some printout
86 ATH_MSG_VERBOSE(" combineEtaPhiPatterns: "
87 <<" eta patterns: " << etaPatternCollection.size()
88 << " phi patterns: " << phiPatternCollection.size()
89 <<std::endl<<"#################################################################################"
90 <<std::endl<<"Print eta pattern collection "<<std::endl<<m_printer->print(etaPatternCollection)
91 <<std::endl<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
92 <<std::endl<<"Print phi pattern collection "<<std::endl<<m_printer->print(phiPatternCollection)
93 );
94
95 for (unsigned int etalevel = 0; etalevel < etaPatternCollection.size(); etalevel++) {
96 CandPrdPatPtr etapattern{etaPatternCollection.at(etalevel), Unowned()};
97 if (etapattern->numberOfContainedPrds() == 0) continue;
98 const Amg::Vector3D etaPatDir = etapattern->globalDirection().unit();
99 const Amg::Vector3D& etaPatPos = etapattern->globalPosition();
100 double theta = etaPatDir.theta();
101 const double phieta = etaPatDir.phi();
102 CxxUtils::sincos scphieta(phieta);
103 CxxUtils::sincos sctheta(theta);
104 double z0 = etaPatPos.z(); // 0 for non -cosmics
105 double rz0 = z0 * sctheta.sn; // closest distance in rz
106
107 double eta_x = etaPatPos.x();
108 double eta_y = etaPatPos.y();
109 double pattern_z0 =z0 + etaPatPos.perp() * sctheta.cs / sctheta.sn;
110 // z belonging to (x0,y0) -> noncosmics just close to
111 // 0 (0.001 * sin/cos)
112 const double eta_r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(eta_x, eta_y, phieta);
113 const double charge = eta_r0>= 0. ? 1. : -1;
114 // Get inverse curvature from eta pattern
115 const double curvature = etapattern->globalDirection().mag();
116 const double invcurvature = curvature > 2 ?charge / curvature : 0.;
117
118 ATH_MSG_DEBUG(" eta pattern info level: " << etalevel << " phi " << phieta << " theta " << theta << " x0 " << eta_x << " y0 "
119 << eta_y << " z0 " << z0 << " hits " << etapattern->numberOfContainedPrds());
120 // flags for cosmics:
121 double min_average_distance = m_maximum_xydistance + m_maximum_rzdistance;
122 CandPrdPatPtr max_phipattern = nullptr;
123 int max_philevel = -1;
124 // flags
125 double dotprodbest = -1.;
126 int phibest = -1;
127 bool ismatched = false;
128 for (unsigned int philevel = 0; philevel < phiPatternCollection.size(); philevel++) {
129 CandPrdPatPtr phipattern{phiPatternCollection.at(philevel), Unowned()};
130 if (phipattern->numberOfContainedPrds() == 0) continue;
131 bool useTightAssociation = false;
132
133 if (m_useTightAssociation && (phipattern->numberOfContainedPrds() > m_maxSizePhiPatternLoose ||
134 etapattern->numberOfContainedPrds() > m_maxSizeEtaPatternLoose)) {
135 if (phipattern->numberOfContainedPrds() > m_maxSizePhiPatternLoose)
136 ATH_MSG_DEBUG(" using tight cuts due to phi hits " << phipattern->numberOfContainedPrds() << " cut "
137 << m_maxSizePhiPatternLoose);
138 if (etapattern->numberOfContainedPrds() > m_maxSizeEtaPatternLoose)
139 ATH_MSG_DEBUG(" using tight cuts due to eta hits " << etapattern->numberOfContainedPrds() << " cut "
140 << m_maxSizeEtaPatternLoose);
141 useTightAssociation = true;
142 }
143
144 const Amg::Vector3D phiPatDir = phipattern->globalDirection().unit();
145 const double dotprod = phiPatDir.dot(etaPatDir);
146
147 if (dotprod > dotprodbest) {
148 dotprodbest = dotprod;
149 phibest = philevel;
150 }
151
152 if (!m_use_cosmics) {
153 ATH_MSG_DEBUG(" eta nr " << etalevel << " phi nr " << philevel << " inproduct " << dotprod << " sin angle "
154 << std::sin(std::acos(dotprod)));
155 }
156
157 double r0{0.}, phipattern_x{0.}, phipattern_y{0.}, phi{phiPatDir.phi()};
158 CxxUtils::sincos scphi{phi};
159 if (m_use_cosmics) {
160 // new information: update parameters
161 std::array<double,4> new_pars = updateParametersForCosmics(*phipattern, *etapattern);
162 r0 = new_pars[0];
163 phi = new_pars[1];
164 rz0 = new_pars[2];
165 theta = new_pars[3];
166 scphi = CxxUtils::sincos(phi);
167 sctheta = CxxUtils::sincos(theta);
168 phipattern_x = r0 * scphi.sn;
169 phipattern_y = r0 * scphi.cs;
170 } else {
171 const Amg::Vector3D posphipattern = phipattern->globalPosition();
172 phipattern_x = posphipattern.x();
173 phipattern_y = posphipattern.y();
174 r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phipattern_x, phipattern_y, phi);
175 }
176
177 // is etapattern on phi pattern?
178 // Loop over hits etapattern
179
180 // int nhits_in_average_eta=0;
181 // int nhits_inside_distance_cut_eta=0;
182
183 if (dotprod <= 0.5 && !m_use_cosmics) continue;
184
185 ATH_MSG_DEBUG(" Matched Eta/phi pattern ");
186
187 // keep track of the number of eta/phi trigger and CSC hits per chamber
188 std::map<Identifier, ChamberInfo> infoPerChamber;
189 std::map<Muon::MuonStationIndex::StIndex, ChamberInfo> infoPerStation;
190 // Loop over hits phi pattern
191 double average_distance{0};
192 int nhits_in_average{0}, nhits_inside_distance_cut{0};
193
194 double phiPatMin{1e9}, phiPatMax{-1e9};
195
196 for (unsigned int phihitid = 0; phihitid < phipattern->numberOfContainedPrds(); phihitid++) {
197 const Trk::PrepRawData* prd = phipattern->prd(phihitid);
198 const Amg::Vector3D& globalposhit = globalPos(prd);
199 double radius_hit = globalposhit.perp();
200 double dotprodradius = sctheta.apply(radius_hit, globalposhit.z());
201 ATH_MSG_VERBOSE("combine hit: " << m_idHelperSvc->toString(prd->identify()) << " dotprod: " << dotprodradius);
202 if (dotprodradius < 0 && !m_use_cosmics) continue;
203
204 double residu_distance_mm{100.*Gaudi::Units::meter};
205 if (m_use_cosmics) {
206 const double perp = scphi.apply(globalposhit.y(), globalposhit.x());
207 residu_distance_mm = MuonHoughMathUtils::signedDistanceToLine(globalposhit.z(), perp, rz0, theta);
208 } else {
209 residu_distance_mm = MuonHoughMathUtils::signedDistanceCurvedToHit(pattern_z0, theta, invcurvature, globalposhit);
210 }
211
212 double distancetoline = std::abs(residu_distance_mm);
213
214 ATH_MSG_VERBOSE(" distance RZ: " << residu_distance_mm);
215 if (distancetoline >= m_maximum_rzdistance) continue;
216
217 ATH_MSG_VERBOSE(" accepted ");
218 nhits_inside_distance_cut++;
219 nhits_in_average++;
220 average_distance += distancetoline;
221
222 if (!useTightAssociation) { continue; }
223 Identifier chId = m_idHelperSvc->chamberId(prd->identify());
224 ChamberInfo& chInfo = infoPerChamber[chId];
225 ++chInfo.nphi;
226 double hitphi = globalposhit.phi();
227 chInfo.phiMin = std::min(hitphi,chInfo.phiMin);
228 chInfo.phiMax = std::min(hitphi,chInfo.phiMax);
229
230 Muon::MuonStationIndex::StIndex stIndex = m_idHelperSvc->stationIndex(prd->identify());
231 ChamberInfo& stInfo = infoPerStation[stIndex];
232 ++stInfo.nphi;
233 stInfo.phiMin = std::min(hitphi, stInfo.phiMin);
234 stInfo.phiMax = std::max(hitphi, stInfo.phiMax);
235 phiPatMin =std::min(hitphi,phiPatMin);
236 phiPatMax =std::max(hitphi,phiPatMax);
237
238 } // size muonpattern
239
240 if (nhits_in_average > 0) average_distance /= nhits_in_average;
241
242 ATH_MSG_DEBUG(" Result for phi pattern: accepted hits " << nhits_inside_distance_cut << " average distance "
243 << average_distance);
244
245 bool etapattern_passed = false;
246 for (unsigned int etahitid = 0; etahitid < etapattern->numberOfContainedPrds(); etahitid++) {
247 const Trk::PrepRawData* prd = etapattern->prd(etahitid);
248 const Amg::Vector3D& etaglobalposhit = globalPos(prd);
249 const double etahitx = etaglobalposhit.x();
250 const double etahity = etaglobalposhit.y();
251
252 if (!m_use_cosmics) {
253 double etadotprod = scphi.apply(etahity, etahitx);
254 // same as in maketruetracksegment
255 if (etadotprod < 0) continue; // hit in wrong hemisphere
256 }
257
258 const double xdiff = phipattern_x - etahitx;
259 const double ydiff = phipattern_y - etahity;
260 const double etahitr = std::hypot(xdiff, ydiff);
261
262 bool hit_passed = false;
263 double etadistancetoline = std::abs(MuonHoughMathUtils::distanceToLine(etahitx, etahity, r0, phi));
264
265 ATH_MSG_VERBOSE("combine: " << m_idHelperSvc->toString(prd->identify()) << " distance xy " << etadistancetoline);
266
267 if (m_use_cosmics) { // phi cone does not work for cosmics since hits
268 // might be close to position of pattern
269
270 const double scale = etahitr / 7000.;
271 hit_passed = etadistancetoline < scale * m_maximum_xydistance;
272 } else if (2 * etadistancetoline < etahitr) { // this corresponds with 30 degrees , normal formula:
273 // (etadistancetoline/etahitr) < sin( Pi/180 * degrees)
274 hit_passed = true;
275 }
276
277 if (!hit_passed) { continue; }
278
279 etapattern_passed = true; // only one hit required
280 // break;
281 ATH_MSG_VERBOSE(" accepted");
282
283 if (!useTightAssociation) { continue ;}
284 Identifier chId = m_idHelperSvc->chamberId(prd->identify());
285
286 ChamberInfo& chInfo = infoPerChamber[chId];
287 ++chInfo.neta;
288 if (m_idHelperSvc->isMdt(prd->identify())) {
289 Muon::MuonStationIndex::StIndex stIndex = m_idHelperSvc->stationIndex(prd->identify());
290 ChamberInfo& stInfo = infoPerStation[stIndex];
291
292 const MuonGM::MdtReadoutElement* mdtDetEl =
293 dynamic_cast<const MuonGM::MdtReadoutElement*>(prd->detectorElement());
294 if (!mdtDetEl) continue;
295
296 const Identifier id = prd->identify();
297 const Trk::Surface& surf = mdtDetEl->surface(id);
298 int layer = m_idHelperSvc->mdtIdHelper().tubeLayer(id);
299 int tube = m_idHelperSvc->mdtIdHelper().tube(id);
300 double halfLength = 0.5 * mdtDetEl->getWireLength(layer, tube);
301 Amg::Vector2D lpLeft(0, -halfLength);
302 Amg::Vector3D gposLeft = surf.localToGlobal(lpLeft);
303 double phiLeft = gposLeft.phi();
304
305 Amg::Vector2D lpRight(0, halfLength);
306 const Amg::Vector3D gposRight = surf.localToGlobal(lpRight);
307 double phiRight = gposRight.phi();
308 double phiMin = std::min(phiRight, phiLeft);
309 double phiMax = std::max(phiLeft ,phiRight);
310
311 Amg::Vector3D tubePos = mdtDetEl->tubePos(id);
312 Amg::Vector3D ROPos = mdtDetEl->ROPos(id);
313 Amg::Vector3D HVPos = 2 * tubePos - ROPos;
314 double tubeL = (HVPos - ROPos).mag();
315 double phiRO = ROPos.phi();
316 double phiHV = HVPos.phi();
317 double phiMinPos = std::min(phiHV, phiRO);
318 double phiMaxPos = std::max(phiRO, phiHV);
319
320 if (std::abs(phiMin - phiMinPos) > 0.01 || std::abs(phiMax - phiMaxPos) > 0.01) {
321 ATH_MSG_DEBUG(" inconsistent Phi!!: from locToGlob (" << phiMin << "," << phiMax << "), from positions ("
322 << phiMinPos << "," << phiMaxPos << ")");
323 }
324 double rotationFraction = 0.;
325 if (phiMin < 0 && phiMax > 0) {
326 if (phiMin < -0.75 * M_PI || phiMax > 0.75 * M_PI)
327 rotationFraction = 1.5;
328 else
329 rotationFraction = 0.5;
330 } else if (phiMax < 0) {
331 rotationFraction = 1.;
332 }
333 double phiMinR = rotatePhi(phiMin, rotationFraction);
334 double phiMaxR = rotatePhi(phiMax, rotationFraction);
335 phiMin = std::min(phiMinR, phiMaxR);
336 phiMax = std::max(phiMinR, phiMaxR);
337
338 phiMinR = rotatePhi(phiMinPos, rotationFraction);
339 phiMaxR = rotatePhi(phiMaxPos, rotationFraction);
340 phiMinPos = std::min(phiMinR, phiMaxR);
341 phiMaxPos = std::max(phiMinR, phiMaxR);
342
343 // enlarge range by 0.1 rad
344 phiMin = phiMin > 0 ? phiMin - 0.1 : phiMin + 0.1;
345 phiMax = phiMax > 0 ? phiMax + 0.1 : phiMax - 0.1;
346
347 double phiMinSec{1.e9},phiMaxSec{-1.e9};
348 if (stInfo.nphi > 0 && stInfo.phiMin < 1000) {
349 phiMinR = rotatePhi(stInfo.phiMin, rotationFraction);
350 phiMaxR = rotatePhi(stInfo.phiMax, rotationFraction);
351 phiMinSec = std::min(phiMinR, phiMaxR);
352 phiMaxSec = std::max(phiMinR, phiMaxR);
353
354 bool inside = true;
355
356 // easy case
357 if (phiMinSec > 0 && phiMaxSec > 0) {
358 if (phiMin > phiMaxSec || phiMax < phiMinSec) inside = false;
359 } else if (phiMinSec < 0 && phiMaxSec < 0) {
360 // easy case (2), always outside
361 inside = false;
362 } else {
363 // finaly check if phiMaxSec large than maxPhi
364 if (phiMax < phiMaxSec) inside = false;
365 }
366 // case with a
367 if (inside) {
368 ++stInfo.ninside;
369 ++chInfo.ninside;
370 ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" Inside ");
371 } else {
372 ++stInfo.noutside;
373 ++chInfo.noutside;
374 ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" Outside ");
375 }
376 }
377
378 phiMinR = rotatePhi(phiPatMin, rotationFraction);
379 phiMaxR = rotatePhi(phiPatMax, rotationFraction);
380 double phiMinPat = std::min(phiMinR, phiMaxR);
381 double phiMaxPat = std::max(phiMinR, phiMaxR);
382
383 bool insidePat = true;
384 // easy case
385 if (phiMinPat > 0 && phiMaxPat > 0) {
386 if (phiMin > phiMaxPat || phiMax < phiMinPat) insidePat = false;
387 } else if (phiMinPat < 0 && phiMaxPat < 0) {
388 // easy case (2), always outside
389 insidePat = false;
390 } else {
391 // finaly check if phiMaxPat large than maxPhi
392 if (phiMax < phiMaxPat) insidePat = false;
393 }
394
395 // case with a
396 if (insidePat) {
397 ++stInfo.ninsidePat;
398 ++chInfo.ninsidePat;
399 ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" InPat ");
400 } else {
401 ++stInfo.noutsidePat;
402 ++chInfo.noutsidePat;
403 ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" OutPat ");
404 }
405 if (myDebug) {
406 ATH_MSG_DEBUG(" : Phi MDT ("
407 << std::setprecision(3) << std::setw(4) << phiMin << "," << std::setw(4) << phiMax << ") "
408 << " from pos (" << std::setprecision(3) << std::setw(4) << phiMinPos << "," << std::setw(4)
409 << phiMaxPos << ") ");
410 if (stInfo.nphi > 0 && stInfo.phiMin < 1000) {
411 ATH_MSG_DEBUG(" phi range (" << std::setprecision(3) << std::setw(4) << stInfo.phiMin << ","
412 << std::setw(4) << stInfo.phiMax << ") ");
413 }
414 ATH_MSG_DEBUG(" pat range (" << std::setprecision(3) << std::setw(4) << phiMinPat << "," << std::setw(4)
415 << phiMaxPat << ") " << m_idHelperSvc->toString(prd->identify()));
416 ATH_MSG_DEBUG(" ATL " << mdtDetEl->getActiveTubeLength(layer, tube) << " WL "
417 << mdtDetEl->getWireLength(layer, tube) << " POSL " << tubeL);
418 }
419 }
420 } // eta pattern
421
422 if (!etapattern_passed) continue; // no etahit close to phipattern, try next phi pattern
423 unsigned int netaPhiPairs = 0;
424 if (useTightAssociation) {
425 // now we have calculated the eta/phi hit content of the 'merged'
426 // pattern
427 for ( auto& [chamberId, chamberInfo] : infoPerChamber) {
428
429 ATH_MSG_VERBOSE(" " << std::setw(32) << m_idHelperSvc->toStringChamber(chamberId) << " eta hits "
430 << chamberInfo.neta << " phi hits " << chamberInfo.nphi << " ninside "
431 << chamberInfo.ninside << " noutside " << chamberInfo.noutside << " ninside "
432 << chamberInfo.ninsidePat << " noutside " << chamberInfo.noutsidePat);
433
434 netaPhiPairs += (chamberInfo.neta && chamberInfo.nphi);
435 }
436
437 ATH_MSG_DEBUG(" eta/phi pattern hit overlap " << netaPhiPairs);
438 if (!etapattern_passed) { ATH_MSG_DEBUG(" failed eta hit match "); }
439 if (nhits_inside_distance_cut < (phipattern->numberOfContainedPrds() * 0.25)) {
440 ATH_MSG_DEBUG(" failed phi hit match ");
441 }
442 if (netaPhiPairs == 0) { ATH_MSG_DEBUG(" Bad match, no overlap "); }
443
444 }
445 ATH_MSG_VERBOSE(" Eta pattern compatible with phi pattern, eta/phi overlap " << netaPhiPairs << " ass phi hits "
446 << nhits_inside_distance_cut << " tot phi hits "
447 << phipattern->numberOfContainedPrds()
448 <<(useTightAssociation ? " using tight association ": "" ));
449
450 // at least 25% matched, to be more robust than 1!
451 if ((!useTightAssociation || netaPhiPairs > 0) && nhits_inside_distance_cut >= (phipattern->numberOfContainedPrds() * 0.25)) {
452 ismatched = true;
453
454 // for non-cosmics try every candidate
455 // for cosmics keep track of best candidate, possible eff loss when phi
456 // patterns are split
457 if (m_use_cosmics) {
458 std::array<double,4> new_pars{r0, phi, pattern_z0, theta};
459 PrdPatternPair updatedpatterns = updatePatternsForCosmics(*phipattern, *etapattern, new_pars);
460 phipattern = std::move(updatedpatterns.first);
461 etapattern = std::move(updatedpatterns.second);
462 ATH_MSG_DEBUG(" Combination accepted with cosmic selection ");
463 } else if (useTightAssociation) {
464 ATH_MSG_DEBUG(" Tight association, cleaning patterns");
465
466 // clean up hits using local phi info
467 std::unique_ptr<Muon::MuonPrdPattern> etaPat = std::make_unique<Muon::MuonPrdPattern>(etapattern->globalPosition(), etapattern->globalDirection());
468 std::unique_ptr<Muon::MuonPrdPattern> phiPat = std::make_unique<Muon::MuonPrdPattern>(phipattern->globalPosition(), phipattern->globalDirection());
469 for (unsigned int etahitid = 0; etahitid < etapattern->numberOfContainedPrds(); ++etahitid) {
470 const Trk::PrepRawData* prd = etapattern->prd(etahitid);
471 const Identifier id = prd->identify();
472 Identifier chId = m_idHelperSvc->chamberId(id);
473 std::map<Identifier, ChamberInfo>::iterator chPos = infoPerChamber.find(chId);
474 if (chPos == infoPerChamber.end()) continue;
475
476 if (m_idHelperSvc->isMdt(id)) {
477 if (chPos->second.ninside == 0 && chPos->second.noutside > 0) continue;
478 if (chPos->second.ninsidePat == 0 && chPos->second.noutsidePat > 0) continue;
479 } else {
480 if (chPos->second.nphi == 0) continue;
481 }
482 etaPat->addPrd(prd);
483 }
484 for (unsigned int phihitid = 0; phihitid < phipattern->numberOfContainedPrds(); ++phihitid) {
485 const Trk::PrepRawData* prd = phipattern->prd(phihitid);
486 const Identifier& id = prd->identify();
487 Identifier chId = m_idHelperSvc->chamberId(id);
488 std::map<Identifier, ChamberInfo>::iterator chPos = infoPerChamber.find(chId);
489 if (chPos == infoPerChamber.end()) continue;
490
491 if (chPos->second.neta == 0) continue;
492 phiPat->addPrd(prd);
493 }
494 phipattern = std::move(phiPat);
495 etapattern = std::move(etaPat);
496 }
497
498 if (!m_bestphimatch) {
499 addCandidate(etapattern, phipattern, candidates, false, phiEtaHitAssMap);
500 ATH_MSG_DEBUG("Candidate FOUND eta " << etalevel << " phi " << philevel << " dotprod: " << dotprod);
501 } else {
502 if (average_distance < min_average_distance) {
503 ATH_MSG_DEBUG(" Selected as best candidate ");
504 min_average_distance = average_distance;
505 max_phipattern = phipattern;
506 max_philevel = philevel;
507 }
508 ATH_MSG_DEBUG(" theta pattern " << etapattern->globalDirection().theta() << " phi "
509 << phipattern->globalDirection().phi() << "average distance " << average_distance
510 << " number of hits " << nhits_inside_distance_cut << " etalevel: " << etalevel);
511 }
512
513 // add recovery for the case we have an inefficiency in the eta hits
514 } else if (useTightAssociation && netaPhiPairs == 0 &&
515 nhits_inside_distance_cut >= (phipattern->numberOfContainedPrds() * 0.25)) {
516 ATH_MSG_DEBUG(" Combination rejected by phi/eta overlap: average distance " << average_distance);
517
518 if (average_distance < min_average_distance) {
519 ATH_MSG_DEBUG(" but selected as best candidate ");
520 min_average_distance = average_distance;
521 max_phipattern = phipattern;
522 max_philevel = philevel;
523 }
524 } // nhits>=25%
525 } // size phi level
526
527 // for cosmics only match best phi pattern
528 if (m_bestphimatch && ismatched) {
529 addCandidate(etapattern, max_phipattern, candidates, true, phiEtaHitAssMap);
530 ATH_MSG_DEBUG("Candidate FOUND eta " << etalevel << " phi " << max_philevel);
531 }
532 // make associated phi pattern for every etapattern and as candidate for
533 // robustness not needed for cosmics when matched, since already done for
534 // first match:
535
536 if (!(m_use_cosmics && m_splitpatterns && ismatched)) {
537 std::unique_ptr<Muon::MuonPrdPattern> assphipattern = makeAssPhiPattern(*etapattern, phiEtaHitAssMap, true);
538 ATH_MSG_DEBUG("No match found, trying to create associated phi pattern ");
539 if (assphipattern) {
540 // make sure ass phi pattern is not a complete subset of one of the
541 // other phi patterns:
542
543 bool subsetcheck = false;
544
545 std::reverse_iterator<std::vector<CandidatePatPair>::iterator> rit =
546 candidates.rbegin();
547 for (; rit != candidates.rend(); ++rit) {
548 if ((*rit).first != etapattern) { break; }
549 // for (;phi_it!=range.second;++phi_it) {
550 if (subset(assphipattern.get(), (*rit).second.get())) {
551 subsetcheck = true;
552 break;
553 }
554 }
555 if (subsetcheck) {
556
557 } else {
558 // these associated phi patterns should be deleted at end of routine:
559
560
561 ATH_MSG_DEBUG("Candidate FOUND eta " << etalevel << " and associated phipattern ");
562
563 ismatched = true;
564
565 // print associated pattern:
566 if (msgLvl(MSG::VERBOSE)) printPattern(assphipattern.get());
567
568 if (m_use_cosmics) {
569 std::array<double,4> new_pars = updateParametersForCosmics(*assphipattern, *etapattern);
570 PrdPatternPair updatedpatterns = updatePatternsForCosmics(*assphipattern, *etapattern, new_pars);
571 assphipattern = std::move(updatedpatterns.first);
572 etapattern = std::move(updatedpatterns.second);
573 }
574 ATH_MSG_DEBUG(" adding eta pattern with recalculated associated phi pattern ");
575
576 addCandidate(etapattern, std::move(assphipattern), candidates, false, phiEtaHitAssMap);
577 }
578 }
579 }
580 if (!ismatched && max_philevel > -1) {
581 addCandidate(etapattern, max_phipattern, candidates, true, phiEtaHitAssMap);
582 ismatched = true;
583 ATH_MSG_DEBUG("No good candidate found, adding best phi pattern " << etalevel << " phi " << max_philevel);
584 }
585 if (!ismatched) {
586 if (msgLvl(MSG::DEBUG)) {
587 ATH_MSG_DEBUG("NO COMBINED Candidate FOUND eta " << etalevel << " phi " << phibest);
588 if (!m_use_cosmics) ATH_MSG_DEBUG("dotprodbest: " << dotprodbest);
589 ATH_MSG_DEBUG("writing out eta pattern (no cleanup)");
590 }
591 candidates.emplace_back(etapattern, nullptr);
592 } else {
593 if (msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG("Candidate was associated to a phi pattern ");
594 }
595 } // size rtheta level
596 return makeCombinedPatterns(candidates);
597
598} // combineEtaPhiPatterns
599
600std::unique_ptr<MuonPrdPatternCollection> MuonCombinePatternTool::makeCombinedPatterns(
601 std::vector<CandidatePatPair>& candidates) const {
602 ATH_MSG_DEBUG("Number of Candidates: " << candidates.size());
603
604 // if (m_use_cosmics == true) {
605 cleanCandidates(candidates);
606 ATH_MSG_DEBUG("Number of Candidates after cleaning: " << candidates.size());
607 //}
608
609 std::unique_ptr<MuonPrdPatternCollection> combinedPatternCollection = std::make_unique<MuonPrdPatternCollection>();
610 int number_comb_patterns = 0;
611 for ( const auto& [etapattern, phipattern] : candidates) {
612 ATH_MSG_DEBUG("Next Candidate");
613 if (phipattern) {
614 std::unique_ptr<Muon::MuonPrdPattern> combinedpattern = makeCombinedPattern(*phipattern, *etapattern);
615 if (combinedpattern) {
616 if (msgLvl(MSG::VERBOSE)) { printPattern(combinedpattern.get()); }
617 combinedPatternCollection->push_back(combinedpattern.release());
618 number_comb_patterns++;
619 } else
620 ATH_MSG_WARNING("combined pattern lost");
621 } else { // no combined match found && no associated phi hits found
622 if (m_splitpatterns && m_use_cosmics) {
623 ATH_MSG_VERBOSE("No combined pattern, eta pattern split based on phi direction of eta pattern ");
624
625 std::vector<PrdPatternPair> splitetapatterns = splitPatternsCylinder(nullptr, etapattern.get());
626 if (splitetapatterns.empty()) {
627 combinedPatternCollection->push_back(etapattern->clone());
628 } else {
629 for (unsigned int i = 0; i < splitetapatterns.size(); i++) {
630 if (splitetapatterns[i].second->numberOfContainedPrds() != 0) {
631 combinedPatternCollection->push_back(splitetapatterns[i].second.release());
632 }
633 }
634 }
635 } else { // don't split pattern
636 ATH_MSG_DEBUG("No combined pattern, eta pattern copied ");
637 combinedPatternCollection->push_back(etapattern->clone());
638 }
639 }
640 }
641
642 ATH_MSG_DEBUG("Number of combined patterns: " << number_comb_patterns << " Number of unmatched etapatterns: "
643 << combinedPatternCollection->size() - number_comb_patterns);
644
645 return combinedPatternCollection;
646}
647
648std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::makeCombinedPattern(const Muon::MuonPrdPattern& phipattern,
649 const Muon::MuonPrdPattern& etapattern) const {
650 ATH_MSG_DEBUG("make combined pattern");
651
652 double phi = phipattern.globalDirection().phi();
653 double theta = etapattern.globalDirection().theta();
654 CxxUtils::sincos scphi(phi);
655 CxxUtils::sincos sctheta(theta);
656
657 double phieta = etapattern.globalDirection().phi();
658 double eta_x = etapattern.globalPosition().x();
659 double eta_y = etapattern.globalPosition().y();
660 double eta_r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(eta_x, eta_y, phieta);
661 double charge = 1.;
662 if (!m_use_cosmics && eta_r0 < 0) charge = -1.;
663 double curvature = etapattern.globalDirection().mag();
664
665 // Store charge in sign of r0 (like in Muon eta pattern)
666
667 const double x0 = charge * (phipattern.globalPosition().x());
668 const double y0 = charge * (phipattern.globalPosition().y());
669 const double z0_phi = m_use_cosmics ?(calculateRz0(etapattern, phi, theta) / sctheta.sn) : 0.; // for non-cosmics
670
671 const Amg::Vector3D dir{curvature * scphi.cs * sctheta.sn, curvature * scphi.sn * sctheta.sn, curvature * sctheta.cs};
672 const Amg::Vector3D pos{x0, y0, z0_phi};
673
674 Muon::MuonPrdPattern::PrdVector comb_prds{};
675 comb_prds.insert(comb_prds.end(),etapattern.prepRawDataVec().begin(),etapattern.prepRawDataVec().end());
676 comb_prds.insert(comb_prds.end(),phipattern.prepRawDataVec().begin(),phipattern.prepRawDataVec().end());
677 std::unique_ptr<Muon::MuonPrdPattern> combinedpattern = std::make_unique<Muon::MuonPrdPattern>(pos, dir, std::move(comb_prds));
678
679 ATH_MSG_DEBUG("Combined pattern with charge " << charge << " curvature " << curvature);
680 ATH_MSG_DEBUG("direction combined pattern: " << scphi.cs * sctheta.sn << " " << scphi.sn * sctheta.sn << " " << sctheta.cs);
681 ATH_MSG_DEBUG("position combined pattern: " << x0 << " " << y0 << " " << z0_phi);
682 ATH_MSG_DEBUG("etapatternsize: " << etapattern.numberOfContainedPrds());
683 ATH_MSG_DEBUG("phipatternsize: " << phipattern.numberOfContainedPrds());
684 ATH_MSG_DEBUG("Combined Track size: " << combinedpattern->numberOfContainedPrds());
685
686 if (m_use_cosmics) {
687 if (msgLvl(MSG::VERBOSE)) ATH_MSG_VERBOSE("No Cleaning for Cosmics!");
688
689 if (msgLvl(MSG::VERBOSE)) { printPattern(combinedpattern.get()); }
690 return combinedpattern;
691 }
692
693 ATH_MSG_DEBUG("Start Cleaning and Recalculating of Combined Pattern");
694
695 bool change = true;
696
697 while (change) {
698 int size_before_cleanup = combinedpattern->numberOfContainedPrds();
699 ATH_MSG_DEBUG("size before cleanup: " << size_before_cleanup);
700
701 std::unique_ptr<Muon::MuonPrdPattern> cleaneduppattern = cleanupCombinedPattern(*combinedpattern);
702
703 int size_after_cleanup = cleaneduppattern->numberOfContainedPrds();
704 ATH_MSG_DEBUG("size after cleanup: " << size_after_cleanup);
705
706 if (size_before_cleanup == size_after_cleanup || size_after_cleanup == 0) {
707 if (msgLvl(MSG::VERBOSE)) { printPattern(cleaneduppattern.get()); }
708 return cleaneduppattern;
709 } else if (size_after_cleanup < size_before_cleanup) {
710 combinedpattern = std::move(cleaneduppattern);
711 } else {
712 change = false;
713 ATH_MSG_FATAL("Cosmic Muon through computer bit? ");
714 }
715 } // while
716 return nullptr;
717} // makeCombinedPattern
718
719std::vector<PrdPatternPair> MuonCombinePatternTool::splitPatterns2D(
720 const Muon::MuonPrdPattern* phipattern, const Muon::MuonPrdPattern* etapattern) const {
721 std::vector<PrdPatternPair> splitPatterns;
722 splitPatterns.reserve(2);
723
724 double phi = phipattern->globalDirection().phi();
725 double theta = etapattern->globalDirection().theta();
726
727 CxxUtils::sincos scphi(phi);
728 CxxUtils::sincos sctheta(theta);
729
730 const Amg::Vector3D dir1{scphi.cs * sctheta.sn, scphi.sn * sctheta.sn, sctheta.cs};
731
732 std::unique_ptr<Muon::MuonPrdPattern> phipattern1 = std::make_unique<Muon::MuonPrdPattern>(phipattern->globalPosition(), dir1); // "lower" pattern (y<0)
733 std::unique_ptr<Muon::MuonPrdPattern> etapattern1 = std::make_unique<Muon::MuonPrdPattern>(etapattern->globalPosition(), dir1); // "lower" pattern (y<0)
734
735 Amg::Vector3D dir2{dir1};
736
737 if (m_flipdirectionforcosmics) {
738 const double newphi = phi + M_PI;
739 const double newtheta = M_PI - etapattern->globalDirection().theta();
740 CxxUtils::sincos scnewphi(newphi);
741 CxxUtils::sincos scnewtheta(newtheta);
742
743 // flip phi and theta for second pattern:
744 dir2 = Amg::Vector3D(scnewphi.cs * scnewtheta.sn, scnewphi.sn * scnewtheta.sn, scnewtheta.cs);
745 }
746
747 std::unique_ptr<Muon::MuonPrdPattern> phipattern2 = std::make_unique<Muon::MuonPrdPattern>(phipattern->globalPosition(), dir2); // "upper" pattern (y>0)
748 std::unique_ptr<Muon::MuonPrdPattern> etapattern2 = std::make_unique<Muon::MuonPrdPattern>(etapattern->globalPosition(), dir2); // "upper" pattern (y>0)
749
750 ATH_MSG_DEBUG(" split pattern1 theta: " << phipattern1->globalDirection().theta() << " phi: " << phipattern1->globalDirection().phi());
751 ATH_MSG_DEBUG(" split pattern2 theta: " << phipattern2->globalDirection().theta() << " phi: " << phipattern2->globalDirection().phi());
752
753 for (unsigned int hitid = 0; hitid < phipattern->numberOfContainedPrds(); hitid++) {
754 const Trk::PrepRawData* prd = phipattern->prd(hitid);
755 const Amg::Vector3D& globalposhit = globalPos(prd);
756 const double dotprod = scphi.apply(globalposhit.y(), globalposhit.x());
757 if (dotprod >= 0) {
758 phipattern1->addPrd(prd);
759 } else {
760 phipattern2->addPrd(prd);
761 }
762 }
763 for (unsigned int hitid = 0; hitid < etapattern->numberOfContainedPrds(); hitid++) {
764 const Trk::PrepRawData* prd = etapattern->prd(hitid);
765 const Amg::Vector3D& globalposhit = globalPos(prd);
766 const double dotprod = scphi.apply(globalposhit.y(), globalposhit.x());
767 if (dotprod >= 0) {
768 etapattern1->addPrd(prd);
769 } else {
770 etapattern2->addPrd(prd);
771 }
772 }
773 //clone: cant have two unique_ptrs to the same object in the container
774 splitPatterns.emplace_back(std::move(phipattern1), std::make_unique<Muon::MuonPrdPattern>(*etapattern1));
775 splitPatterns.emplace_back(std::move(etapattern1), std::move(etapattern2));
776 return splitPatterns;
777}
778
779std::vector<PrdPatternPair> MuonCombinePatternTool::splitPatterns3D(
780 const Muon::MuonPrdPattern* phipattern, const Muon::MuonPrdPattern* etapattern) const {
781 // phi pattern may be 0
782
783 std::vector<PrdPatternPair> splitPatterns;
784 splitPatterns.reserve(2);
785
786 // phi and etapattern should have identical directions and positions (except
787 // for z_0), but new built anyway to be robust
788
789 Amg::Vector3D globalpos{Amg::Vector3D::Zero()};
790 if (phipattern) {
791 globalpos = phipattern->globalPosition();
792 } else {
793 globalpos = etapattern->globalPosition();
794 }
795
796 double phi;
797 if (phipattern) {
798 phi = phipattern->globalDirection().phi();
799 } else {
800 phi = etapattern->globalDirection().phi();
801 }
802
803 const double theta = etapattern->globalDirection().theta();
804
805 CxxUtils::sincos scphi(phi);
806 CxxUtils::sincos sctheta(theta);
807
808 Amg::Vector3D dir1{scphi.cs * sctheta.sn, scphi.sn * sctheta.sn, sctheta.cs};
809
810 Amg::Vector3D dir2 = dir1;
811 if (m_flipdirectionforcosmics) {
812 const double newphi = phi + M_PI;
813 const double newtheta = M_PI - theta;
814
815 CxxUtils::sincos scnewphi(newphi);
816 CxxUtils::sincos scnewtheta(newtheta);
817
818 // flip phi and theta for second pattern:
819 dir2 = Amg::Vector3D(scnewphi.cs * scnewtheta.sn, scnewphi.sn * scnewtheta.sn, scnewtheta.cs);
820 }
821
822 std::unique_ptr<Muon::MuonPrdPattern> etapattern1 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir1); //
823 std::unique_ptr<Muon::MuonPrdPattern> etapattern2 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir2); // "upper" pattern (y>0)
824 std::unique_ptr<Muon::MuonPrdPattern> phipattern1, phipattern2{};
825
826 if (phipattern) {
827 phipattern1 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir1); // "lower" pattern (y<0)
828 phipattern2 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir2); // "upper" pattern (y>0)
829 }
830
831
832
833 if (msgLvl(MSG::DEBUG)) {
834 ATH_MSG_DEBUG(" split pattern theta: " << theta << " phi: " << phi);
835 ATH_MSG_DEBUG(" split pattern1 theta: " << etapattern1->globalDirection().theta()
836 << " phi: " << etapattern1->globalDirection().phi());
837 ATH_MSG_DEBUG(" split pattern2 theta: " << etapattern2->globalDirection().theta()
838 << " phi: " << etapattern2->globalDirection().phi());
839 Amg::Vector3D splitpoint = MuonHoughMathUtils::shortestPointOfLineToOrigin(globalpos, phi, theta);
840 ATH_MSG_DEBUG(" splitpoint, x: " << splitpoint[0] << " y: " << splitpoint[1] << " z: " << splitpoint[2]);
841 }
842
843 double d_x = scphi.cs * sctheta.sn;
844 double d_y = scphi.sn * sctheta.sn;
845
846 if (phipattern) {
847 for (unsigned int hitid = 0; hitid < phipattern->numberOfContainedPrds(); hitid++) {
848 const Trk::PrepRawData* prd = phipattern->prd(hitid);
849 const Amg::Vector3D& globalposhit = globalPos(prd);
850 const double hitx = globalposhit.x();
851 const double hity = globalposhit.y();
852 const double hitz = globalposhit.z();
853 const double dotprod = hitx * d_x + hity * d_y + hitz * sctheta.cs;
854 if (dotprod >= 0) {
855 phipattern1->addPrd(prd);
856 } else {
857 phipattern2->addPrd(prd);
858 }
859 ATH_MSG_VERBOSE(" dotprod: " << dotprod);
860 }
861 }
862
863 for (unsigned int hitid = 0; hitid < etapattern->numberOfContainedPrds(); hitid++) {
864 const Trk::PrepRawData* prd = etapattern->prd(hitid);
865 const Amg::Vector3D& globalposhit = globalPos(prd);
866 const double hitx = globalposhit.x();
867 const double hity = globalposhit.y();
868 const double hitz = globalposhit.z();
869 const double dotprod = hitx * d_x + hity * d_y + hitz * sctheta.cs;
870 if (dotprod >= 0) {
871 etapattern1->addPrd(prd);
872 } else {
873 etapattern2->addPrd(prd);
874 }
875 ATH_MSG_VERBOSE(" dotprod: " << dotprod);
876 }
877
878 if (msgLvl(MSG::DEBUG)) {
879 if (phipattern) {
880 ATH_MSG_DEBUG("Final size, phi: " << phipattern1->numberOfContainedPrds() << " " << phipattern2->numberOfContainedPrds());
881 }
882 ATH_MSG_DEBUG("Final size, eta: " << etapattern1->numberOfContainedPrds() << " " << etapattern2->numberOfContainedPrds());
883 }
884 splitPatterns.emplace_back(std::move(phipattern1), std::move(etapattern1));
885 splitPatterns.emplace_back(std::move(phipattern2), std::move(etapattern2));
886
887 return splitPatterns;
888}
889
890std::vector<PrdPatternPair> MuonCombinePatternTool::splitPatternsCylinder(
891 const Muon::MuonPrdPattern* phipattern, const Muon::MuonPrdPattern* etapattern) const {
892 // phi pattern may be 0
893 Amg::Vector3D patternpos{Amg::Vector3D::Zero()};
894 double phi{0.};
895
896 if (phipattern) {
897 patternpos = phipattern->globalPosition();
898 phi = phipattern->globalDirection().phi();
899 } else {
900 patternpos = etapattern->globalPosition();
901 phi = etapattern->globalDirection().phi();
902 }
903
904 const double theta = etapattern->globalDirection().theta();
905
906 // decide if track is split:
907
908 // calculate intersection with pattern and calorimeter cylinder (r0=4000,
909 // c0=6000) if there is one, then track will be split
910
911 bool intersect = MuonHoughMathUtils::lineThroughCylinder(patternpos, phi, theta, MuonHough::radius_cylinder, MuonHough::z_cylinder);
912
913 if (!intersect) { // no split
914 ATH_MSG_DEBUG("Pattern not through calorimeter -> do not split ");
915 return {};
916 }
917
918 return splitPatterns3D(phipattern, etapattern);
919}
920
921std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::cleanupCombinedPattern(const Muon::MuonPrdPattern& combinedpattern) const {
922 const double phipattern = combinedpattern.globalDirection().phi();
923 const double thetapattern = combinedpattern.globalDirection().theta();
924
925 CxxUtils::sincos scthetapattern(thetapattern);
926
927 const Amg::Vector3D& patternpos = combinedpattern.globalPosition();
928 const double posx = patternpos.x();
929 const double posy = patternpos.y();
930 const double posz = patternpos.z();
931
932 double invcurvature = 0.;
933 double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(posx, posy, phipattern);
934 double charge = 1.;
935 if (r0 < 0) charge = -1.;
936 double curvature = combinedpattern.globalDirection().mag();
937 if (curvature > 2) invcurvature = charge / curvature;
938
939 ATH_MSG_DEBUG("cleaned up pattern: phi " << phipattern << " theta: " << thetapattern << " position: " << posx << " " << posy << " "
940 << posz);
941 ATH_MSG_DEBUG("Cleanup pattern charge " << charge << " curvature " << curvature);
942
943 std::unique_ptr<Muon::MuonPrdPattern> combinedpattern_cleaned =
944 std::make_unique<Muon::MuonPrdPattern>(combinedpattern.globalPosition(), combinedpattern.globalDirection());
945
946 for (unsigned int hitid = 0; hitid < combinedpattern.numberOfContainedPrds(); hitid++) {
947 const Trk::PrepRawData* prd = combinedpattern.prd(hitid);
948 const Amg::Vector3D& globalposhit = globalPos(prd);
949
950 double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(posx, posy, phipattern);
951 double distance_xy = MuonHoughMathUtils::distanceToLine(globalposhit.x(), globalposhit.y(), r0, phipattern);
952
953 double radius_pattern = globalposhit.perp();
954 double z0 = posz - radius_pattern * scthetapattern.cs / scthetapattern.sn;
955
956 double distance_rz = MuonHoughMathUtils::signedDistanceCurvedToHit(z0, thetapattern, invcurvature, globalposhit);
957
958 const double scale = std::max(1., globalposhit.mag() / 7000.);
959 ATH_MSG_VERBOSE("hit: " << globalposhit);
960 ATH_MSG_VERBOSE("CLEAN distancetopattern: "
961 << " dist xy " << distance_xy << " dist rz " << distance_rz << " scale: " << scale);
962 if (std::abs(distance_xy) < scale * m_maximum_xydistance && std::abs(distance_rz) < m_maximum_rzdistance) {
963 combinedpattern_cleaned->addPrd(prd);
964 } else {
965 ATH_MSG_DEBUG("Hit discarded: " << hitid << " dis xy " << distance_xy << " dis rz " << distance_rz);
966 ATH_MSG_DEBUG("Hit info: "<<m_idHelperSvc->toString(prd->identify()));
967 }
968 }
969
970 ATH_MSG_DEBUG("size of cleaned pattern: " << combinedpattern_cleaned->numberOfContainedPrds());
971
972 if (combinedpattern_cleaned->numberOfContainedPrds() == 0 && combinedpattern.numberOfContainedPrds() != 0) {
973 ATH_MSG_DEBUG(
974 "cleaned up pattern is empty (should happen only when initially no phi "
975 "pattern found and phi hits are added by ascociation map)");
976 }
977
978 return combinedpattern_cleaned;
979}
980
981std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::makeAssPhiPattern(
982 const Muon::MuonPrdPattern& muonpattern,
984 const EtaPhiHitAssocMap& phiEtaHitAssMap,
985 bool check_already_on_pattern) const {
986 // bool hits_added = false;
987 const unsigned int size = muonpattern.numberOfContainedPrds();
988
989 PrepDataSet hits;
990 if (check_already_on_pattern) {
991 for (unsigned int i = 0; i < size; i++) { hits.insert(muonpattern.prd(i)); }
992 }
993 std::vector<const Trk::PrepRawData*> phihits;
994 for (unsigned int i = 0; i < size; i++) {
995 const Trk::PrepRawData* prd = muonpattern.prd(i);
996 // check on type of prd?
997 const Muon::MuonCluster* muoncluster = dynamic_cast<const Muon::MuonCluster*>(prd);
998 if (!muoncluster) continue;
999 EtaPhiHitAssocMap::const_iterator itr = phiEtaHitAssMap.find(prd);
1001 if (itr == phiEtaHitAssMap.end()) {
1002 continue;
1003 }
1004 std::copy_if(itr->second.begin(), itr->second.end(), std::back_inserter(phihits),
1005 [&check_already_on_pattern, &hits](const Trk::PrepRawData* phiHit){
1006 return !check_already_on_pattern || hits.insert(phiHit).second;
1007 });
1008 }
1009
1010
1011 if (phihits.empty()) { return nullptr; }
1012
1013 std::unique_ptr<Muon::MuonPrdPattern> phipattern;
1014
1015 double phi = 0., sin_phi = 0., cos_phi = 0.;
1016 if (m_use_cosmics) {
1017 phi = muonpattern.globalDirection().phi();
1018 } else {
1019 for (const Trk::PrepRawData* phiHit : phihits) {
1020 const Amg::Vector3D& globalposhit = globalPos(phiHit);
1021 CxxUtils::sincos scphihit(globalposhit.phi());
1022 sin_phi += scphihit.sn;
1023 cos_phi += scphihit.cs;
1024 }
1025 phi = std::atan2(sin_phi, cos_phi);
1026 }
1027
1028 const double curvature = muonpattern.globalDirection().mag();
1029 const double theta = muonpattern.globalDirection().theta();
1030 CxxUtils::sincos scphi(phi);
1031 CxxUtils::sincos sctheta(theta);
1032
1033 if (m_use_cosmics) {
1034 phipattern = std::make_unique<Muon::MuonPrdPattern>(muonpattern.globalPosition(), muonpattern.globalDirection(), phihits);
1035 } else {
1036 static const Amg::Vector3D globalpos{0.001, 0.001, 0.};
1037 const Amg::Vector3D& globaldir{curvature * scphi.cs * sctheta.sn,
1038 curvature * scphi.sn * sctheta.sn, curvature * sctheta.cs};
1039 phipattern = std::make_unique<Muon::MuonPrdPattern>(globalpos, globaldir, phihits);
1040 }
1041
1042
1043 // perform cleaning on newly created phipattern:
1044 std::unique_ptr<Muon::MuonPrdPattern> phipatternclean = cleanPhiPattern(std::move(phipattern));
1045
1046 if (phipatternclean->numberOfContainedPrds() <= 0) {
1047 phipatternclean.reset();
1048 }
1049 return phipatternclean;
1050}
1051
1052std::array<double,4> MuonCombinePatternTool::updateParametersForCosmics(const Muon::MuonPrdPattern& phipattern,
1053 const Muon::MuonPrdPattern& etapattern) const {
1054 // method returns r0, phi, rz0, theta
1055
1056 const unsigned int etasize = etapattern.numberOfContainedPrds();
1057 const unsigned int phisize = phipattern.numberOfContainedPrds();
1058
1059 const Amg::Vector3D& etaglobalpos = etapattern.globalPosition();
1060 const Amg::Vector3D& etaglobaldir = etapattern.globalDirection();
1061
1062 const Amg::Vector3D& phiglobalpos = phipattern.globalPosition();
1063 const Amg::Vector3D& phiglobaldir = phipattern.globalDirection();
1064
1065 std::array<double,4> old_pars{0};
1066
1067 old_pars[0] = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phiglobaldir.phi());
1068 old_pars[1] = phiglobaldir.phi();
1069
1070 const double theta_orig = etaglobaldir.theta();
1071 old_pars[2] = etaglobalpos.z() * std::sin(theta_orig);
1072 old_pars[3] = theta_orig;
1073
1074 if (phisize + etasize <= 1) return old_pars;
1075
1076 // first calculate new phi, r0 (initial value , -Pi/2):
1077 std::pair<double, double> rphi_start = calculateR0Phi(phipattern, etapattern);
1078 // for stabilising (can be cpu-optimised greatly):
1079 std::pair<double, double> rphi = calculateR0Phi(phipattern, etapattern, rphi_start.first);
1080
1081 if (msgLvl(MSG::DEBUG) && std::abs(std::sin(rphi.first - rphi_start.first)) > 0.15 &&
1082 std::abs(std::sin(etaglobaldir.phi() - phiglobaldir.phi())) < 0.15) {
1083 ATH_MSG_DEBUG("unexpected phi change!");
1084 ATH_MSG_DEBUG("phi first: " << rphi_start.first << " phi second: " << rphi.first);
1085 ATH_MSG_DEBUG("phi etapattern: " << etaglobaldir.phi() << " phi phipattern: " << phiglobaldir.phi());
1086 }
1087
1088 const double phi = rphi.first;
1089 const double r0 = rphi.second;
1090
1091 CxxUtils::sincos scphi(phi);
1092
1093 // calculate new theta and rz0: (not using phi hits this time, as eta pattern
1094 // is larger in general
1095
1096 double av_radii{0.}, av_z{0.};
1097
1098 for (unsigned int i = 0; i < etasize; ++i) {
1099 const Trk::PrepRawData* prd = etapattern.prd(i);
1100 const Amg::Vector3D globalposhit = globalPos(prd);
1101 av_radii += scphi.apply(globalposhit.y(), globalposhit.x());
1102 av_z += globalposhit.z();
1103 }
1104
1105 if (etasize > 0) {
1106 av_radii /= etasize;
1107 av_z /= etasize;
1108 }
1109 double sumr = 0.;
1110 double sumz = 0.;
1111 for (unsigned int i = 0; i < etasize; i++) {
1112 const Trk::PrepRawData* prd = etapattern.prd(i);
1113 const Amg::Vector3D& globalposhit = globalPos(prd);
1114
1115 double radius = scphi.apply(globalposhit.y(), globalposhit.x()); // hitx*scphi.cs + hity*scphi.sn;
1116 double r_offset = radius - av_radii;
1117 double z_offset = globalposhit.z() - av_z;
1118 double weight = r_offset * r_offset + z_offset * z_offset;
1119 int sign = 1;
1120 if (r_offset * std::cos(theta_orig) + z_offset * std::sin(theta_orig) < 0) { sign = -1; }
1121 sumr += weight * sign * r_offset;
1122 sumz += weight * sign * z_offset;
1123 }
1124
1125 // const double sum_tan = sum_tanweight/sum_weight;
1126
1127 ATH_MSG_VERBOSE("av_z : " << av_z << " av_radii: " << av_radii << " sumr: " << sumr << " sumz: " << sumz);
1128 if (std::abs(sumr) < 0.000001 || std::abs(sumz) < 0.000001) return old_pars;
1129
1130 double theta = std::atan2(sumr, sumz);
1131
1132 if (theta < 0) theta += M_PI;
1133
1134 double rz0 = calculateRz0(etapattern, phi, theta);
1135
1136 // ATH_MSG_DEBUG("old method rz0: " << sctheta.apply(av_z,-av_radii) );
1137 // const double rz0 = sctheta.apply(av_z,-av_radii); // (av_z * sctheta.sn) -
1138 // av_radii * sctheta.cs;
1139
1140 std::array<double,4> new_pars {r0,phi,rz0,theta};
1141
1142 ATH_MSG_VERBOSE("updated parameters: r0: " << new_pars[0] << " phi: " << new_pars[1] << " rz0: " << new_pars[2]
1143 << " theta: " << new_pars[3]);
1144 ATH_MSG_VERBOSE("old parameters: r0: " << old_pars[0] << " phi: " << old_pars[1] << " rz0: " << old_pars[2] << " theta: " << old_pars[3]);
1145
1146 if (msgLvl(MSG::VERBOSE)) {
1147 ATH_MSG_VERBOSE("phisize: " << phisize << " etasize: " << etasize);
1148 for (unsigned int i = 0; i < phisize; i++) {
1149 const Trk::PrepRawData* prd = phipattern.prd(i);
1150 const Amg::Vector3D& globalposhit = globalPos(prd);
1151 double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), r0, phi);
1152 ATH_MSG_VERBOSE("distance to updated parameters in xy: " << distance);
1153 distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), old_pars[0], old_pars[1]);
1154 ATH_MSG_VERBOSE("old distance phi hit: " << distance);
1155 }
1156 for (unsigned int i = 0; i < etasize; i++) {
1157 const Trk::PrepRawData* prd = etapattern.prd(i);
1158 const Amg::Vector3D& globalposhit = globalPos(prd);
1159 double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), r0, phi);
1160 ATH_MSG_VERBOSE("distance to updated parameters in xy: " << distance);
1161 distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), old_pars[0], old_pars[1]);
1162 ATH_MSG_VERBOSE("old distance eta hit: " << distance);
1163 }
1164 for (unsigned int i = 0; i < etasize; ++i) {
1165 const Trk::PrepRawData* prd = etapattern.prd(i);
1166 const Amg::Vector3D& globalposhit = globalPos(prd);
1167 double perp = scphi.apply(globalposhit.y(),
1168 globalposhit.x()); // globalposhit.x()*scphi.cs + globalposhit.y()*scphi.sn;
1169 double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit[Amg::z], perp, rz0, theta);
1170 ATH_MSG_VERBOSE("distance to updated parameters in Rz: " << distance);
1171 distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit[Amg::z], perp, old_pars[2], old_pars[3]);
1172 ATH_MSG_VERBOSE("old distance: " << distance);
1173 }
1174 }
1175 return new_pars;
1176}
1177
1178std::pair<double, double> MuonCombinePatternTool::calculateR0Phi(const Muon::MuonPrdPattern& phipattern,
1179 const Muon::MuonPrdPattern& etapattern, double phi_est) const {
1180 // use eta pattern as well, since phi patterns consist sometimes of only 1
1181 // station etahit error 200mrad (2Pi/16*2), phi hit 20mrad (should be hough
1182 // binsize (18 mrad), but prefer ~ factor 10 for stabilility)
1183
1184 // test if lever_arm > 2 m before updating , if not old values are used
1185
1186 ATH_MSG_VERBOSE("calculateR0Phi");
1187
1188 CxxUtils::sincos scphi_est(phi_est);
1189
1190 const unsigned int etasize = etapattern.numberOfContainedPrds();
1191 const unsigned int phisize = phipattern.numberOfContainedPrds();
1192
1193 const Amg::Vector3D& etaglobaldir = etapattern.globalDirection();
1194 const double phi_etapattern = etaglobaldir.phi();
1195 CxxUtils::sincos scphi_eta(phi_etapattern);
1196
1197 const Amg::Vector3D& phiglobaldir = phipattern.globalDirection();
1198 const Amg::Vector3D& phiglobalpos = phipattern.globalPosition();
1199 const double phi_phipattern = phiglobaldir.phi();
1200 CxxUtils::sincos scphi_phi(phi_phipattern);
1201
1202 const double phi_error_inv = 1. / 20.;
1203 const double phi_error_inv2 = phi_error_inv * phi_error_inv;
1204 const double eta_error_inv = 1. / 400.;
1205 const double eta_error_inv2 = eta_error_inv * eta_error_inv;
1206
1207 // from MuonHoughPattern::updateParametersRPhi (partial code duplication.. :(
1208 // )
1209
1210 double sum_etax{0.}, sum_etay{0.}, sum_phix{0.}, sum_phiy{0.};
1211
1212 // calculate average point
1213
1214 for (unsigned int i = 0; i < etasize; i++) {
1215 const Trk::PrepRawData* prd = etapattern.prd(i);
1216 const Amg::Vector3D& globalposhit = globalPos(prd);
1217 sum_etax += globalposhit.x();
1218 sum_etay += globalposhit.y();
1219 }
1220
1221 for (unsigned int i = 0; i < phisize; i++) {
1222 const Trk::PrepRawData* prd = phipattern.prd(i);
1223 const Amg::Vector3D& globalposhit = globalPos(prd);
1224 sum_phix += globalposhit.x();
1225 sum_phiy += globalposhit.y();
1226 }
1227
1228 const double av_x = (eta_error_inv2 * sum_etax + phi_error_inv2 * sum_phix) / (eta_error_inv2 * etasize + phi_error_inv2 * phisize);
1229 const double av_y = (eta_error_inv2 * sum_etay + phi_error_inv2 * sum_phiy) / (eta_error_inv2 * etasize + phi_error_inv2 * phisize);
1230
1231 ATH_MSG_VERBOSE(" av_x: " << av_x << " av_y: " << av_y);
1232
1233 // calculate weighted sum:
1234
1235 double sumx {0.}, sumy {0.};
1236
1237 // keep track of extreme points
1238
1239 double x_min {0.}, x_max {0.}, y_min {0.}, y_max {0.}, lever_min {0.}, lever_max{0.};
1240
1241 for (unsigned int i = 0; i < etasize; i++) {
1242 const Trk::PrepRawData* prd = etapattern.prd(i);
1243 const Amg::Vector3D& globalposhit = globalPos(prd);
1244 double x_offset = globalposhit.x() - av_x;
1245 double y_offset = globalposhit.y() - av_y;
1246 double height_squared = x_offset * x_offset + y_offset * y_offset;
1247 double weight = height_squared * eta_error_inv2;
1248 int sign = 1;
1249 if (x_offset * scphi_est.cs + y_offset * scphi_est.sn < 0) { sign = -1; }
1250 sumx += weight * sign * x_offset;
1251 sumy += weight * sign * y_offset;
1252
1253 if (sign == 1 && height_squared > lever_max) {
1254 lever_max = height_squared;
1255 x_max = globalposhit.x();
1256 y_max = globalposhit.y();
1257 } else if (sign == -1 && height_squared > lever_min) {
1258 lever_min = height_squared;
1259 x_min = globalposhit.x();
1260 y_min = globalposhit.y();
1261 }
1262 }
1263
1264 for (unsigned int i = 0; i < phisize; i++) {
1265 const Trk::PrepRawData* prd = phipattern.prd(i);
1266 const Amg::Vector3D& globalposhit = globalPos(prd);
1267 double x_offset = globalposhit.x() - av_x;
1268 double y_offset = globalposhit.y() - av_y;
1269 double height_squared = x_offset * x_offset + y_offset * y_offset;
1270 double weight = height_squared * phi_error_inv2;
1271 int sign = 1;
1272 if (x_offset * scphi_est.cs + y_offset * scphi_est.sn < 0) { sign = -1; }
1273 sumx += weight * sign * x_offset;
1274 sumy += weight * sign * y_offset;
1275
1276 if (sign == 1 && height_squared > lever_max) {
1277 lever_max = height_squared;
1278 x_max = globalposhit.x();
1279 y_max = globalposhit.y();
1280 } else if (sign == -1 && height_squared > lever_min) {
1281 lever_min = height_squared;
1282 x_min = globalposhit.x();
1283 y_min = globalposhit.y();
1284 }
1285 }
1286
1287 ATH_MSG_VERBOSE("av_x : " << av_x << " av_y: " << av_y << " sumx: " << sumx << " sumy: " << sumy);
1288
1289 if (std::abs(sumx) < 0.000001 || std::abs(sumy) < 0.000001) {
1290 ATH_MSG_DEBUG(" sum too small to update");
1291
1292 return std::make_pair(phi_phipattern, MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phi_phipattern));
1293 }
1294
1295 // lever arm has to be larger than 2 m, else no update:
1296 if (std::hypot(x_max - x_min , y_max - y_min) < 2000) {
1297 ATH_MSG_VERBOSE("lever arm too small: av_x : " << std::sqrt((x_max - x_min) * (x_max - x_min) + (y_max - y_min) * (y_max - y_min))
1298 << " x_max: " << x_max << " x_min: " << x_min << " y_max: " << y_max
1299 << " y_min: " << y_min);
1300 return std::make_pair(phi_phipattern, MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phi_phipattern));
1301 }
1302
1303 double phi_fit = std::atan2(sumy, sumx);
1304 if (phi_fit > 0) phi_fit -= M_PI; // phi between 0,-Pi for cosmics!
1305 CxxUtils::sincos scphi(phi_fit);
1306 const double r0_fit = scphi.apply(av_x, -av_y); // av_x * scphi.sn - av_y * scphi.cs;
1307
1308 return std::make_pair(phi_fit, r0_fit);
1309}
1310
1311double MuonCombinePatternTool::calculateRz0(const Muon::MuonPrdPattern& pattern, double phi, double theta) {
1312 double nhits = pattern.numberOfContainedPrds();
1313 CxxUtils::sincos sctheta(theta);
1314 CxxUtils::sincos scphi(phi);
1315
1316 /*
1317 x = r_0 sin(phi) + t*cos(phi)sin(theta)
1318 y = - r_0 cos(phi) + t*sin(phi)sin(theta)
1319 z = z_0 + t*cos(theta)
1320
1321 2 methods (average point of average_z0):
1322
1323 r^2 = x^2+y^2 = r_0^2 + t^2*sin^2(theta)
1324 (this projects the hit radially onto the line)
1325
1326 Not so good as the radius of the hits can be smaller than r_0
1327
1328 method based on:
1329
1330 x*cos(phi) + y*sin(phi) = t*sin(theta)
1331 (projects the hit on the line in x-y and calculates the distance to r_0)
1332
1333 works best
1334 */
1335
1336 // method 3:
1337
1338 double rz0 = 0.;
1339 for (unsigned int i = 0; i < nhits; i++) {
1340 const Trk::PrepRawData* prd = pattern.prd(i);
1341 const Amg::Vector3D& poshit = globalPos(prd);
1342 int sign = (poshit.x() * scphi.cs + poshit.y() * scphi.sn > 0) ? 1 : -1;
1343 rz0 += poshit.z() * sctheta.sn - sign * sctheta.cs * poshit.perp();
1344 }
1345
1346 if (nhits > 0) rz0 /= nhits;
1347 return rz0;
1348}
1349
1350MuonCombinePatternTool::PrdPatternPair MuonCombinePatternTool::updatePatternsForCosmics(
1351 const Muon::MuonPrdPattern& phipattern, const Muon::MuonPrdPattern& etapattern, const std::array<double,4>& new_pars) const {
1352 double phi = new_pars[1];
1353 double theta = new_pars[3];
1354
1355 CxxUtils::sincos scphi(phi);
1356 CxxUtils::sincos sctheta(theta);
1357
1358 double x0 = new_pars[0] * scphi.sn;
1359 double y0 = -new_pars[0] * scphi.cs;
1360 double z0_phi = new_pars[2];
1361 double z0_eta = new_pars[2];
1362 if (std::abs(sctheta.sn) > 1e-7) {
1363 z0_phi = (new_pars[2] + new_pars[0] * sctheta.cs) / sctheta.sn; // z0 belonging to (x0,y0)
1364 z0_eta = new_pars[2] / sctheta.sn; // z0 of rz0
1365 }
1366
1367 const Amg::Vector3D globalDir{scphi.cs * sctheta.sn, scphi.sn * sctheta.sn, sctheta.cs};
1368 const Amg::Vector3D globalPosPhi{x0, y0, z0_phi};
1369 const Amg::Vector3D globalPosEta{x0, y0, z0_eta};
1370 ATH_MSG_VERBOSE("updatePatternsForCosmics() -- eta pattern "<<std::endl
1371 <<m_printer->print(*etapattern.prd(0))<<std::endl
1372 <<" -- phi pattern "<<std::endl
1373 <<m_printer->print(*phipattern.prd(0))<<std::endl
1374 <<"Theta: " <<theta<<" globalPosEta: "<<Amg::toString(globalPosEta)<<", globalDir: "<<Amg::toString(globalDir));
1375
1376
1377 std::unique_ptr<Muon::MuonPrdPattern> updatedphipattern = std::make_unique<Muon::MuonPrdPattern>(globalPosPhi, globalDir, phipattern.prepRawDataVec());
1378 std::unique_ptr<Muon::MuonPrdPattern> updatedetapattern = std::make_unique<Muon::MuonPrdPattern>(globalPosEta, globalDir, etapattern.prepRawDataVec());
1379
1380 return std::make_pair(std::move(updatedphipattern), std::move(updatedetapattern));
1381}
1382
1383std::unique_ptr<MuonPatternCombinationCollection> MuonCombinePatternTool::makePatternCombinations(const MuonPrdPatternCollection& muonpatterns) const {
1384 ATH_MSG_DEBUG("makePatternCombinations");
1385
1386 std::unique_ptr<MuonPatternCombinationCollection> patterncombinations = std::make_unique<MuonPatternCombinationCollection>();
1387
1388 for (const Muon::MuonPattern* pit : muonpatterns) {
1389 const Amg::Vector3D roadmom = pit->globalDirection();
1390 const Amg::Vector3D roadpos = pit->globalPosition();
1391 ATH_MSG_DEBUG("phi: " << roadmom.phi() << " eta: " << roadmom.eta());
1392 ATH_MSG_DEBUG("x: " << roadpos.x() << " y: " << roadpos.y() << " z: " << roadpos.z());
1393
1394 // sort pattern per chamber
1395 std::map<Identifier, std::vector<const Trk::PrepRawData*>> chamberMap;
1396 for (unsigned int i = 0; i < pit->numberOfContainedPrds(); ++i) {
1397 const Trk::PrepRawData* prd = pit->prd(i);
1398 Identifier channelId = prd->identify();
1399 const Identifier moduleId = m_idHelperSvc->chamberId(channelId);
1400 std::vector<const Trk::PrepRawData*>& chambVec = chamberMap[moduleId];
1401 if (chambVec.empty()) chambVec.reserve(pit->numberOfContainedPrds());
1402 chambVec.push_back(prd);
1403 }
1404
1405 // build chamberintersect vector
1406 std::vector<Muon::MuonPatternChamberIntersect> mpciVec;
1407 mpciVec.reserve(chamberMap.size());
1408 Amg::Vector3D patpose{Amg::Vector3D::Zero()}, patdire{Amg::Vector3D::Zero()};
1409 for ( const auto& [moduleId, chambVec] : chamberMap) {
1410 const Trk::PrepRawData* prd = chambVec.front();
1411 const Amg::Vector3D& globalpos = globalPos(prd);
1412 if (m_use_cosmics) {
1413 // not flip
1414 patdire = roadmom.unit();
1415 patpose = roadpos;
1416 } else {
1417 MuonHoughMathUtils::extrapolateCurvedRoad(roadpos, roadmom, globalpos, patpose, patdire);
1418 }
1419
1420 Muon::MuonPatternChamberIntersect mpci = Muon::MuonPatternChamberIntersect(patpose, patdire, chambVec);
1421 mpciVec.push_back(mpci);
1422 }
1423
1424 ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" Pattern position "<<Amg::toString(roadpos));
1425 Amg::Vector3D tvertex{Amg::Vector3D::Zero()};
1426 Trk::TrackParameters* parameters = new Trk::Perigee(roadpos, roadmom, 1., tvertex); // if -1 then charge flipped anyway
1427 Muon::MuonPatternCombination* combination = new Muon::MuonPatternCombination(parameters, mpciVec);
1428 patterncombinations->push_back(combination);
1429 }
1430 return patterncombinations;
1431}
1432
1433bool MuonCombinePatternTool::subset(const Muon::MuonPrdPattern* pattern1, const Muon::MuonPrdPattern* pattern2) {
1435
1436 // first check if pattern 1 is not larger than 2:
1437 if (pattern1->numberOfContainedPrds() > pattern2->numberOfContainedPrds()) { return false; }
1438
1439 PrepDataSet hits;
1440 for (unsigned int hitnr = 0; hitnr < pattern2->numberOfContainedPrds(); ++hitnr) { hits.insert(pattern2->prd(hitnr)); }
1441
1442 for (unsigned int hitnr = 0; hitnr < pattern1->numberOfContainedPrds(); ++hitnr) {
1443 if (!hits.count(pattern1->prd(hitnr))) {
1444 return false;
1445 }
1446 }
1447 return true;
1448}
1449
1450bool MuonCombinePatternTool::subset(std::pair<PrepDataSet,
1451 PrepDataSet>& candidate1,
1452 std::pair<PrepDataSet,
1453 PrepDataSet>& candidate2) {
1455
1456 // first check if pattern 1 is not larger than 2:
1457 if (candidate1.first.size() > candidate2.first.size() ||
1458 candidate1.second.size() > candidate2.second.size()) { return false; }
1459 for (const Trk::PrepRawData* find_me : candidate1.first) {
1460 if (!candidate2.first.count(find_me)) { return false; }
1461 }
1462 for (const Trk::PrepRawData* find_me : candidate1.second) {
1463 if (candidate2.second.count(find_me)) { return false; }
1464 }
1465
1466 return true;
1467}
1468
1469void MuonCombinePatternTool::printPattern(const Muon::MuonPrdPattern* muonpattern) const {
1470 if (msgLvl(MSG::VERBOSE)) {
1471 ATH_MSG_VERBOSE("Printout of Pattern: ");
1472
1473 const Amg::Vector3D& pos = muonpattern->globalPosition();
1474 const Amg::Vector3D& dir = muonpattern->globalDirection();
1475
1476 ATH_MSG_VERBOSE("pos: x: " << pos.x() << " y: " << pos.y() << " z: " << pos.z());
1477 ATH_MSG_VERBOSE("dir: x: " << dir.x() << " y: " << dir.y() << " z: " << dir.z());
1478 ATH_MSG_VERBOSE("phi: " << dir.phi() << " theta: " << dir.theta() << " rz0: " << pos.z() * std::sin(dir.theta()));
1479
1480 for (unsigned int k = 0; k < muonpattern->numberOfContainedPrds(); k++) {
1481 const Trk::PrepRawData* prd = muonpattern->prd(k);
1482 const Muon::MdtPrepData* mdtprd = dynamic_cast<const Muon::MdtPrepData*>(prd);
1483 if (mdtprd) {
1484 const Trk::Surface& surface = mdtprd->detectorElement()->surface(mdtprd->identify());
1485 const Amg::Vector3D& gpos = surface.center();
1486 ATH_MSG_VERBOSE("mdt " << k << " x: " << gpos.x() << " y: " << gpos.y() << " z: " << gpos.z());
1487 } else if (!mdtprd) {
1488 const Muon::MuonCluster* muoncluster = dynamic_cast<const Muon::MuonCluster*>(prd);
1489 if (muoncluster) {
1490 const Amg::Vector3D& gpos = muoncluster->globalPosition();
1491 ATH_MSG_VERBOSE("cluster " << k << " x: " << gpos.x() << " y: " << gpos.y() << " z: " << gpos.z());
1492 }
1493 if (!muoncluster) { ATH_MSG_VERBOSE("no muon prd? "); }
1494 }
1495 }
1496 }
1497}
1498
1499std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::cleanPhiPattern(std::unique_ptr<Muon::MuonPrdPattern> phipattern) const {
1500 const Amg::Vector3D& olddir = phipattern->globalDirection();
1501 const double theta = olddir.theta();
1502 const unsigned int size = phipattern->numberOfContainedPrds();
1503
1504 if (msgLvl(MSG::DEBUG)) {
1505 ATH_MSG_DEBUG("Start Phi hits cleaning with " << size << " hits "
1506 << " theta " << theta);
1507 const Amg::Vector3D& oldpos = phipattern->globalPosition();
1508 double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(oldpos.x(), oldpos.y(), olddir.phi());
1509 ATH_MSG_DEBUG("Start Phi: " << olddir.phi() << " r0: " << r0);
1510 }
1511
1512 // need internal class to be able to remove hits fast
1513 std::unique_ptr<MuonHoughPattern> newpattern = std::make_unique<MuonHoughPattern>(MuonHough::hough_xy);
1514 for (unsigned int phihitnr = 0; phihitnr < size; phihitnr++) { newpattern->addHit(std::make_shared<MuonHoughHit>(phipattern->prd(phihitnr))); }
1515 newpattern->updateParametersRPhi(m_use_cosmics);
1516 double phi = newpattern->getEPhi();
1517 double r0 = newpattern->getERPhi();
1518
1519 CxxUtils::sincos scphi(phi);
1520
1521 constexpr int number_of_iterations = 4;
1522 std::array<double,number_of_iterations> cutvalues {1000., 500., 250., 125.};
1523 if (m_use_cosmics) {
1524 cutvalues[0] = 5000.;
1525 cutvalues[1] = 2500.;
1526 cutvalues[2] = 1250.;
1527 cutvalues[3] = 1250.;
1528 }
1529
1530 for (int it = 0; it < number_of_iterations; ++it) {
1531 ATH_MSG_VERBOSE("iteration " << it << " cutvalue: " << cutvalues[it]);
1532 bool change = true;
1533 while (change) {
1534 ATH_MSG_VERBOSE("size: " << newpattern->size() << " r0: " << r0 << " phi: " << phi);
1535
1536 double max_dist = 0.;
1537 unsigned int max_i = 99999;
1538 for (unsigned int i = 0; i < newpattern->size(); i++) {
1539 double dist =
1540 scphi.apply(newpattern->getHitx(i), -newpattern->getHity(i)) - r0; // newpattern->getHitx(i) * scphi.sn -
1541 // newpattern->getHity(i) * scphi.cs - r0;
1542 ATH_MSG_VERBOSE("Dist: " << dist);
1543 if (std::fabs(dist) > std::abs(max_dist)) {
1544 max_dist = dist;
1545 max_i = i;
1546 }
1547 }
1548 if (std::abs(max_dist) < cutvalues[it]) {
1549 change = false;
1550 } else {
1551 newpattern->removeHit(max_i);
1552 newpattern->updateParametersRPhi(m_use_cosmics);
1553 phi = newpattern->getEPhi();
1554 r0 = newpattern->getERPhi();
1555 scphi = CxxUtils::sincos(phi);
1556 }
1557 }
1558 }
1559
1560 ATH_MSG_DEBUG("Final size: " << newpattern->size() << " r0: " << r0 << " phi: " << phi);
1561
1562 // update parameters rz (not very important as later overwritten)
1563 // put r0 to IP for collisions
1564
1565 double thetanew = 0.;
1566 double r0_new = 1.; // put 1 mm r0 value
1567
1568 if (m_use_cosmics) { r0_new = r0; }
1569
1570 unsigned int nPatterns = newpattern->size();
1571 for (unsigned int i = 0; i < nPatterns; i++) { thetanew += newpattern->getTheta(i); }
1572
1573 if (nPatterns > 0) thetanew /= nPatterns;
1574
1575 double z0_new = 0.;
1576
1577 double x0_new = r0_new * scphi.sn;
1578 double y0_new = -r0_new * scphi.cs;
1579 CxxUtils::sincos sctheta(thetanew);
1580
1581 const Amg::Vector3D pos {x0_new, y0_new, z0_new};
1582 const Amg::Vector3D dir {sctheta.sn * scphi.cs, sctheta.sn * scphi.sn, sctheta.cs};
1583
1584 std::unique_ptr<Muon::MuonPrdPattern> cleanpattern = std::make_unique<Muon::MuonPrdPattern>(pos, dir);
1585
1586 for (unsigned int i = 0; i < newpattern->size(); i++) { cleanpattern->addPrd(newpattern->getPrd(i)); }
1587
1588 return cleanpattern;
1589}
1590
1591void MuonCombinePatternTool::addCandidate(const CandPrdPatPtr& etapattern, const CandPrdPatPtr& phipattern,
1592 std::vector<CandidatePatPair>& candidates, bool add_asspattern,
1593 const EtaPhiHitAssocMap& phiEtaHitAssMap) const {
1594 if (!m_use_cosmics || !m_splitpatterns) {
1595 candidates.emplace_back(etapattern, phipattern);
1596 return;
1597 }
1598
1599 std::vector<PrdPatternPair> splitpatterns = splitPatternsCylinder(phipattern.get(), etapattern.get());
1600
1601 if (splitpatterns.empty()) {
1602 candidates.emplace_back(etapattern, phipattern);
1603 }
1604
1605 else {
1606 for (auto& [phiPattern, etaPattern] : splitpatterns) {
1607 // skip when empty eta pattern , possible duplication when associated phi
1608 // pattern is found, but then will be cleaned later
1609 if (etaPattern->numberOfContainedPrds() == 0) {
1610 continue;
1611 }
1612 candidates.emplace_back(std::move(etaPattern), std::move(phiPattern));
1613 }
1614 }
1615
1616 // make associated pattern don't split eta pattern yet, but split based on phi
1617 // of ass. pattern bool asspattern_added = false;
1618 if (!add_asspattern) { return ;}
1619 std::unique_ptr<Muon::MuonPrdPattern> assphipattern = makeAssPhiPattern(*etapattern, phiEtaHitAssMap, true);
1620
1621 if (!assphipattern) { return;}
1622
1623 // print associated pattern:
1624 if (msgLvl(MSG::VERBOSE)) {
1625 ATH_MSG_VERBOSE("Associated Pattern: ");
1626 printPattern(assphipattern.get());
1627 }
1628
1629 // update parameters:
1630 std::array<double,4> new_pars = updateParametersForCosmics(*assphipattern, *etapattern);
1631 PrdPatternPair updatedpatterns = updatePatternsForCosmics(*assphipattern, *etapattern, new_pars);
1632
1633 std::unique_ptr<Muon::MuonPrdPattern>& cosmicPhiPattern = updatedpatterns.first;
1634 std::unique_ptr<Muon::MuonPrdPattern>& cosmicEtaPattern = updatedpatterns.second;
1635
1636 std::vector<PrdPatternPair> splitpatterns_ass = splitPatternsCylinder(cosmicPhiPattern.get(), cosmicEtaPattern.get());
1637 if (splitpatterns_ass.empty()) {
1638 candidates.emplace_back(std::move(cosmicEtaPattern), std::move(cosmicPhiPattern));
1639 return;
1640 }
1641
1642
1643 for (auto& [splitPhiPattern, splitEtaPattern] : splitpatterns_ass) {
1644 if (splitPhiPattern->numberOfContainedPrds() == 0 ||
1645 splitEtaPattern->numberOfContainedPrds() == 0) {
1646 continue;
1647 }
1648 candidates.emplace_back(std::move(splitEtaPattern), std::move(splitPhiPattern));
1649 }
1650}
1651
1652void MuonCombinePatternTool::cleanCandidates(std::vector<CandidatePatPair>& candidates) {
1653 std::vector<CandidatePatPair>::iterator it1;
1654 std::vector<CandidatePatPair>::iterator it2;
1655 // map between set of prd's (eta and phi) and candidates , stored for speed
1656 std::map<CandidatePatPair, std::pair<PrepDataSet, PrepDataSet>> hitsMap;
1657
1658 // fill map
1659 for (it1 = candidates.begin(); it1 != candidates.end(); ++it1) {
1660 PrepDataSet etahits;
1661 for (unsigned int hitnr = 0; hitnr < (*it1).first->numberOfContainedPrds(); hitnr++) { etahits.insert((*it1).first->prd(hitnr)); }
1662 PrepDataSet phihits;
1663 if ((*it1).second) { // phi pattern might be 0!
1664 for (unsigned int hitnr = 0; hitnr < (*it1).second->numberOfContainedPrds(); hitnr++) {
1665 phihits.insert((*it1).second->prd(hitnr));
1666 }
1667 }
1668 hitsMap.insert(std::make_pair((*it1), std::make_pair(etahits, phihits)));
1669 }
1670 // cppcheck-suppress invalidContainer; candidate.erase(it2) does not invalidate it1
1671 for (it1 = candidates.begin(); it1 != candidates.end(); ++it1) {
1672 std::pair<PrepDataSet, PrepDataSet>& hits1 = hitsMap[(*it1)];
1673 it2 = it1 + 1;
1674 while (it2 != candidates.end()) {
1675 std::pair<PrepDataSet, PrepDataSet>& hits2 = hitsMap[(*it2)];
1676 if (subset((hits2), (hits1))) { // 2 subset of 1, remove 2 // in case of
1677 // equality best (earliest) is kept
1678 it2 = candidates.erase(it2); // it2 points to next item, it1 not invalidated!
1679 } else if (subset((hits1), (hits2))) { // 1 subset of 2, remove 1
1680 it1 = candidates.erase(it1); // it1 points to next item, it2 invalidated!
1681 it2 = it1 + 1;
1682 //cppcheck-suppress selfAssignment
1683 hits1 = hitsMap[(*it1)]; // redefine hits1
1684 } else {
1685 ++it2;
1686 }
1687 }
1688 }
1689}
M_PI
#define M_PI
Definition ActiveFraction.h:11
perp
Scalar perp() const
perp method - perpendicular length
Definition AmgMatrixBasePlugin.h:44
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
mag
Scalar mag() const
mag method
Definition AmgMatrixBasePlugin.h:26
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
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
charge
double charge(const T &p)
Definition AtlasPID.h:997
PrdPatternPair
MuonCombinePatternTool::PrdPatternPair PrdPatternPair
Definition MuonCombinePatternTool.cxx:36
MuonPrdPatternCollection
DataVector< Muon::MuonPrdPattern > MuonPrdPatternCollection
Definition MuonPatternCollection.h:20
sign
int sign(int a)
Definition TRT_StrawNeighbourSvc.h:108
y
#define y
x
#define x
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition AthCommonMsg.h:30
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
MuonCombinePatternTool::printPattern
void printPattern(const Muon::MuonPrdPattern *muonpattern) const
print out pattern hits
Definition MuonCombinePatternTool.cxx:1469
MuonCombinePatternTool::m_maximum_xydistance
const double m_maximum_xydistance
distance cut in xy for hits
Definition MuonCombinePatternTool.h:135
MuonCombinePatternTool::addCandidate
void addCandidate(const CandPrdPatPtr &etapattern, const CandPrdPatPtr &phipattern, std::vector< CandidatePatPair > &candidates, bool add_asspattern, const EtaPhiHitAssocMap &phiEtaHitAssMap) const
adds eta,phi pair to candidate vector, also performs splitting and associated pattern (only for cosmi...
Definition MuonCombinePatternTool.cxx:1591
MuonCombinePatternTool::makePatternCombinations
virtual std::unique_ptr< MuonPatternCombinationCollection > makePatternCombinations(const MuonPrdPatternCollection &muonpatterns) const override
converts MuonPrdPatterns into MuonPatternCombinationCollection MuonPatternCombinationCollection are d...
Definition MuonCombinePatternTool.cxx:1383
MuonCombinePatternTool::m_printer
PublicToolHandle< Muon::MuonEDMPrinterTool > m_printer
Definition MuonCombinePatternTool.h:163
MuonCombinePatternTool::cleanCandidates
static void cleanCandidates(std::vector< CandidatePatPair > &candidates)
clean candidates from subsets or duplicates
Definition MuonCombinePatternTool.cxx:1652
MuonCombinePatternTool::makeCombinedPatterns
std::unique_ptr< MuonPrdPatternCollection > makeCombinedPatterns(std::vector< CandidatePatPair > &candidates) const
make combined pattern from all candidates, removes duplicates with phi when no overlap with eta patte...
Definition MuonCombinePatternTool.cxx:600
MuonCombinePatternTool::m_use_cosmics
bool m_use_cosmics
use cosmic settings
Definition MuonCombinePatternTool.h:140
MuonCombinePatternTool::CandPrdPatPtr
std::shared_ptr< const Muon::MuonPrdPattern > CandPrdPatPtr
Definition MuonCombinePatternTool.h:35
MuonCombinePatternTool::makeAssPhiPattern
std::unique_ptr< Muon::MuonPrdPattern > makeAssPhiPattern(const Muon::MuonPrdPattern &pattern, const EtaPhiHitAssocMap &phiEtaHitAssMap, bool check=false) const
make combined phi pattern by associating phi hits to noncombined eta pattern, return 0 if no phi meas...
Definition MuonCombinePatternTool.cxx:981
MuonCombinePatternTool::m_nodiscarding
bool m_nodiscarding
don't discard any candidates based on gasgap assocation (true by default)
Definition MuonCombinePatternTool.h:147
MuonCombinePatternTool::m_bestphimatch
bool m_bestphimatch
take only best phi match as candidate or take all phi matches (false by default, but true for cosmics...
Definition MuonCombinePatternTool.h:151
MuonCombinePatternTool::m_splitpatterns
bool m_splitpatterns
split patterns (only for cosmics)
Definition MuonCombinePatternTool.h:143
MuonCombinePatternTool::updatePatternsForCosmics
PrdPatternPair updatePatternsForCosmics(const Muon::MuonPrdPattern &phipattern, const Muon::MuonPrdPattern &etapattern, const std::array< double, 4 > &new_pars) const
update patterns based on new track parameters (used only for cosmics) builds 2 new prd patterns
Definition MuonCombinePatternTool.cxx:1350
MuonCombinePatternTool::calculateR0Phi
std::pair< double, double > calculateR0Phi(const Muon::MuonPrdPattern &phipattern, const Muon::MuonPrdPattern &etapattern, double phi_estimate=-M_PI_2) const
calculate phi and r0 for cosmic patterns, phi estimate needs to be given
Definition MuonCombinePatternTool.cxx:1178
MuonCombinePatternTool::calculateRz0
static double calculateRz0(const Muon::MuonPrdPattern &pattern, double phi, double theta)
calculate rz0 for cosmic pattern
Definition MuonCombinePatternTool.cxx:1311
MuonCombinePatternTool::PrdPatternPair
std::pair< std::unique_ptr< Muon::MuonPrdPattern >, std::unique_ptr< Muon::MuonPrdPattern > > PrdPatternPair
Definition MuonCombinePatternTool.h:33
MuonCombinePatternTool::subset
static bool subset(const Muon::MuonPrdPattern *pattern1, const Muon::MuonPrdPattern *pattern2)
is pattern1 a complete subset of other pattern2?
Definition MuonCombinePatternTool.cxx:1433
MuonCombinePatternTool::m_flipdirectionforcosmics
bool m_flipdirectionforcosmics
flip direction for cosmics after splitting as if coming from IP (false by default)
Definition MuonCombinePatternTool.h:155
MuonCombinePatternTool::cleanupCombinedPattern
std::unique_ptr< Muon::MuonPrdPattern > cleanupCombinedPattern(const Muon::MuonPrdPattern &combinedpattern) const
clean combined pattern, remove outliers
Definition MuonCombinePatternTool.cxx:921
MuonCombinePatternTool::m_maximum_rzdistance
const double m_maximum_rzdistance
distance cut in rz for hits
Definition MuonCombinePatternTool.h:137
MuonCombinePatternTool::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition MuonCombinePatternTool.h:161
MuonCombinePatternTool::MuonCombinePatternTool
MuonCombinePatternTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition MuonCombinePatternTool.cxx:38
MuonCombinePatternTool::combineEtaPhiPatterns
std::unique_ptr< MuonPrdPatternCollection > combineEtaPhiPatterns(const MuonPrdPatternCollection &phiPatternCollection, const MuonPrdPatternCollection &etaPatternCollection, const EtaPhiHitAssocMap &phiEtaHitAssMap) const override
Combines phi and eta pattern collection into a new combined pattern collection.
Definition MuonCombinePatternTool.cxx:69
MuonCombinePatternTool::initialize
StatusCode initialize() override
Definition MuonCombinePatternTool.cxx:58
MuonCombinePatternTool::splitPatterns2D
std::vector< PrdPatternPair > splitPatterns2D(const Muon::MuonPrdPattern *phipattern, const Muon::MuonPrdPattern *etapattern) const
split patterns in two at point closest to IP in rphi
Definition MuonCombinePatternTool.cxx:719
MuonCombinePatternTool::makeCombinedPattern
std::unique_ptr< Muon::MuonPrdPattern > makeCombinedPattern(const Muon::MuonPrdPattern &phipattern, const Muon::MuonPrdPattern &etapattern) const override
Combines phi and eta pattern into a new combined pattern.
Definition MuonCombinePatternTool.cxx:648
MuonCombinePatternTool::splitPatternsCylinder
std::vector< PrdPatternPair > splitPatternsCylinder(const Muon::MuonPrdPattern *phipattern, const Muon::MuonPrdPattern *etapattern) const
split patterns in two when crossing calorimeter at point closest to IP in 3D (should be same as split...
Definition MuonCombinePatternTool.cxx:890
MuonCombinePatternTool::splitPatterns3D
std::vector< PrdPatternPair > splitPatterns3D(const Muon::MuonPrdPattern *phipattern, const Muon::MuonPrdPattern *etapattern) const
split patterns in two at point closest to IP in 3D
Definition MuonCombinePatternTool.cxx:779
MuonCombinePatternTool::updateParametersForCosmics
std::array< double, 4 > updateParametersForCosmics(const Muon::MuonPrdPattern &phipattern, const Muon::MuonPrdPattern &etapattern) const
calculate new track parameters of match (only for cosmics!) returns [r0, phi, rz0,...
Definition MuonCombinePatternTool.cxx:1052
MuonCombinePatternTool::cleanPhiPattern
std::unique_ptr< Muon::MuonPrdPattern > cleanPhiPattern(std::unique_ptr< Muon::MuonPrdPattern > phipattern) const
clean phi pattern, similar as in MuonHoughPatternTool, used for newly created phi patterns based on h...
Definition MuonCombinePatternTool.cxx:1499
MuonGM::MdtReadoutElement::getActiveTubeLength
double getActiveTubeLength(const int tubeLayer, const int tube) const
MuonGM::MdtReadoutElement::tubePos
Amg::Vector3D tubePos(const Identifier &id) const
Returns the global position of the given tube.
MuonGM::MdtReadoutElement::ROPos
Amg::Vector3D ROPos(const int tubelayer, const int tube) const
Definition MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:288
MuonGM::MdtReadoutElement::surface
virtual const Trk::Surface & surface() const override final
Return surface associated with this detector element.
Definition MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:885
MuonGM::MdtReadoutElement::getWireLength
double getWireLength(const int tubeLayer, const int tube) const
MuonHoughMathUtils::lineThroughCylinder
static bool lineThroughCylinder(const Amg::Vector3D &vec, double phi, double theta, double r_0, double z_0)
calculates if line (x,y,z,phi,theta) crosses cylinder (r_0,z_0) around origin
Definition MuonHoughMathUtils.cxx:161
MuonHoughMathUtils::signedDistanceOfLineToOrigin2D
static double signedDistanceOfLineToOrigin2D(double x, double y, double phi)
signed distance of line with point (x,y) and angle phi to origin
Definition MuonHoughMathUtils.cxx:110
MuonHoughMathUtils::shortestPointOfLineToOrigin
static Amg::Vector3D shortestPointOfLineToOrigin(const Amg::Vector3D &vec, double phi, double theta)
calculates the 3d-point closest to origin
Definition MuonHoughMathUtils.cxx:139
MuonHoughMathUtils::signedDistanceCurvedToHit
static double signedDistanceCurvedToHit(double z0, double theta, double invcurvature, const Amg::Vector3D &hit)
calculates distance of point (x,y,z) to curved track with z0, theta and invcurvature for curved track...
Definition MuonHoughMathUtils.cxx:207
MuonHoughMathUtils::distanceToLine
static double distanceToLine(double x0, double y0, double r0, double phi)
distance from (x0,y0) to the line (r0,phi), phi in rad
Definition MuonHoughMathUtils.cxx:40
MuonHoughMathUtils::signedDistanceToLine
static double signedDistanceToLine(double x0, double y0, double r0, double phi)
distance from (x0,y0) to the line (r0,phi), phi in rad
Definition MuonHoughMathUtils.cxx:31
MuonHoughMathUtils::extrapolateCurvedRoad
static void extrapolateCurvedRoad(const Amg::Vector3D &roadpos, const Amg::Vector3D &roadmom, const Amg::Vector3D &pos, Amg::Vector3D &roadpose, Amg::Vector3D &roaddire)
extrapolates road to global position
Definition MuonHoughMathUtils.cxx:264
Muon::IMuonCombinePatternTool::PrepDataSet
std::set< const Trk::PrepRawData *, IdentifierPrdLess > PrepDataSet
Definition IMuonCombinePatternTool.h:24
Muon::IMuonCombinePatternTool::EtaPhiHitAssocMap
std::map< const Trk::PrepRawData *, PrepDataSet > EtaPhiHitAssocMap
Definition IMuonCombinePatternTool.h:25
Muon::MdtPrepData
Class to represent measurements from the Monitored Drift Tubes.
Definition MdtPrepData.h:33
Muon::MdtPrepData::detectorElement
virtual const MuonGM::MdtReadoutElement * detectorElement() const override
Returns the detector element corresponding to this PRD.
Definition MdtPrepData.h:141
Muon::MuonCluster::globalPosition
virtual const Amg::Vector3D & globalPosition() const =0
Returns the global position of the measurement (calculated on the fly)
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::MuonPattern
Basic class for patterns in the muon spectrometer consistig out of a list of Trk::PrepRawData objects...
Definition MuonPattern.h:23
Muon::MuonPattern::globalDirection
const Amg::Vector3D & globalDirection() const
Global direction of the pattern.
Definition MuonPattern.h:58
Muon::MuonPattern::globalPosition
const Amg::Vector3D & globalPosition() const
Global position of the pattern.
Definition MuonPattern.h:57
Muon::MuonPrdPattern
Class to store a pattern in the muon system containing PrepRawData pointers.
Definition MuonPrdPattern.h:27
Muon::MuonPrdPattern::numberOfContainedPrds
virtual unsigned int numberOfContainedPrds() const
Number or PrepRawData contained by this Pattern.
Definition MuonPrdPattern.h:64
Muon::MuonPrdPattern::prd
virtual const Trk::PrepRawData * prd(unsigned int index) const
returns the PrepRawData objects depending on the integer, return zero if index out of range
Definition MuonPrdPattern.h:66
Muon::MuonPrdPattern::prepRawDataVec
const PrdVector & prepRawDataVec() const
Definition MuonPrdPattern.h:73
Muon::MuonPrdPattern::PrdVector
std::vector< const Trk::PrepRawData * > PrdVector
Definition MuonPrdPattern.h:29
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::identify
Identifier identify() const
return the identifier
Trk::PrepRawData::type
virtual bool type(PrepRawDataType type) const
Interface method checking the type.
Definition PrepRawData.h:133
Trk::Surface
Abstract Base Class for tracking surfaces.
Definition Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
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.
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition GeoPrimitivesToStringConverter.h:40
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition GeoPrimitives.h:48
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
MuonHough::radius_cylinder
constexpr double radius_cylinder
radius of cylinder
Definition MuonHoughMathUtils.h:24
MuonHough::z_cylinder
constexpr double z_cylinder
length of cylinder
Definition MuonHoughMathUtils.h:26
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition MuonStationIndex.h:23
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
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
CxxUtils::sincos::apply
double apply(double a, double b) const
Definition sincos.h:57
type
Generated on for ATLAS Offline Software by doxygen 1.14.0