ATLAS Offline Software
MuonCombinePatternTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 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"
17 namespace{
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 }
36 using PrdPatternPair = MuonCombinePatternTool::PrdPatternPair;
37 
38 MuonCombinePatternTool::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 
58 StatusCode 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 
69 std::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 
600 std::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 
648 std::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 
719 std::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  splitPatterns.emplace_back(std::move(phipattern1), std::move(etapattern1));
774  splitPatterns.emplace_back(std::move(etapattern1), std::move(etapattern2));
775  return splitPatterns;
776 }
777 
778 std::vector<PrdPatternPair> MuonCombinePatternTool::splitPatterns3D(
779  const Muon::MuonPrdPattern* phipattern, const Muon::MuonPrdPattern* etapattern) const {
780  // phi pattern may be 0
781 
782  std::vector<PrdPatternPair> splitPatterns;
783  splitPatterns.reserve(2);
784 
785  // phi and etapattern should have identical directions and positions (except
786  // for z_0), but new built anyway to be robust
787 
788  Amg::Vector3D globalpos{Amg::Vector3D::Zero()};
789  if (phipattern) {
790  globalpos = phipattern->globalPosition();
791  } else {
792  globalpos = etapattern->globalPosition();
793  }
794 
795  double phi;
796  if (phipattern) {
797  phi = phipattern->globalDirection().phi();
798  } else {
799  phi = etapattern->globalDirection().phi();
800  }
801 
802  const double theta = etapattern->globalDirection().theta();
803 
804  CxxUtils::sincos scphi(phi);
805  CxxUtils::sincos sctheta(theta);
806 
807  Amg::Vector3D dir1{scphi.cs * sctheta.sn, scphi.sn * sctheta.sn, sctheta.cs};
808 
809  Amg::Vector3D dir2 = dir1;
810  if (m_flipdirectionforcosmics) {
811  const double newphi = phi + M_PI;
812  const double newtheta = M_PI - theta;
813 
814  CxxUtils::sincos scnewphi(newphi);
815  CxxUtils::sincos scnewtheta(newtheta);
816 
817  // flip phi and theta for second pattern:
818  dir2 = Amg::Vector3D(scnewphi.cs * scnewtheta.sn, scnewphi.sn * scnewtheta.sn, scnewtheta.cs);
819  }
820 
821  std::unique_ptr<Muon::MuonPrdPattern> etapattern1 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir1); //
822  std::unique_ptr<Muon::MuonPrdPattern> etapattern2 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir2); // "upper" pattern (y>0)
823  std::unique_ptr<Muon::MuonPrdPattern> phipattern1, phipattern2{};
824 
825  if (phipattern) {
826  phipattern1 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir1); // "lower" pattern (y<0)
827  phipattern2 = std::make_unique<Muon::MuonPrdPattern>(globalpos, dir2); // "upper" pattern (y>0)
828  }
829 
830 
831 
832  if (msgLvl(MSG::DEBUG)) {
833  ATH_MSG_DEBUG(" split pattern theta: " << theta << " phi: " << phi);
834  ATH_MSG_DEBUG(" split pattern1 theta: " << etapattern1->globalDirection().theta()
835  << " phi: " << etapattern1->globalDirection().phi());
836  ATH_MSG_DEBUG(" split pattern2 theta: " << etapattern2->globalDirection().theta()
837  << " phi: " << etapattern2->globalDirection().phi());
838  Amg::Vector3D splitpoint = MuonHoughMathUtils::shortestPointOfLineToOrigin(globalpos, phi, theta);
839  ATH_MSG_DEBUG(" splitpoint, x: " << splitpoint[0] << " y: " << splitpoint[1] << " z: " << splitpoint[2]);
840  }
841 
842  double d_x = scphi.cs * sctheta.sn;
843  double d_y = scphi.sn * sctheta.sn;
844 
845  if (phipattern) {
846  for (unsigned int hitid = 0; hitid < phipattern->numberOfContainedPrds(); hitid++) {
847  const Trk::PrepRawData* prd = phipattern->prd(hitid);
848  const Amg::Vector3D& globalposhit = globalPos(prd);
849  const double hitx = globalposhit.x();
850  const double hity = globalposhit.y();
851  const double hitz = globalposhit.z();
852  const double dotprod = hitx * d_x + hity * d_y + hitz * sctheta.cs;
853  if (dotprod >= 0) {
854  phipattern1->addPrd(prd);
855  } else {
856  phipattern2->addPrd(prd);
857  }
858  ATH_MSG_VERBOSE(" dotprod: " << dotprod);
859  }
860  }
861 
862  for (unsigned int hitid = 0; hitid < etapattern->numberOfContainedPrds(); hitid++) {
863  const Trk::PrepRawData* prd = etapattern->prd(hitid);
864  const Amg::Vector3D& globalposhit = globalPos(prd);
865  const double hitx = globalposhit.x();
866  const double hity = globalposhit.y();
867  const double hitz = globalposhit.z();
868  const double dotprod = hitx * d_x + hity * d_y + hitz * sctheta.cs;
869  if (dotprod >= 0) {
870  etapattern1->addPrd(prd);
871  } else {
872  etapattern2->addPrd(prd);
873  }
874  ATH_MSG_VERBOSE(" dotprod: " << dotprod);
875  }
876 
877  if (msgLvl(MSG::DEBUG)) {
878  if (phipattern) {
879  ATH_MSG_DEBUG("Final size, phi: " << phipattern1->numberOfContainedPrds() << " " << phipattern2->numberOfContainedPrds());
880  }
881  ATH_MSG_DEBUG("Final size, eta: " << etapattern1->numberOfContainedPrds() << " " << etapattern2->numberOfContainedPrds());
882  }
883  splitPatterns.emplace_back(std::move(phipattern1), std::move(etapattern1));
884  splitPatterns.emplace_back(std::move(phipattern2), std::move(etapattern2));
885 
886  return splitPatterns;
887 }
888 
889 std::vector<PrdPatternPair> MuonCombinePatternTool::splitPatternsCylinder(
890  const Muon::MuonPrdPattern* phipattern, const Muon::MuonPrdPattern* etapattern) const {
891  // phi pattern may be 0
892  Amg::Vector3D patternpos{Amg::Vector3D::Zero()};
893  double phi{0.};
894 
895  if (phipattern) {
896  patternpos = phipattern->globalPosition();
897  phi = phipattern->globalDirection().phi();
898  } else {
899  patternpos = etapattern->globalPosition();
900  phi = etapattern->globalDirection().phi();
901  }
902 
903  const double theta = etapattern->globalDirection().theta();
904 
905  // decide if track is split:
906 
907  // calculate intersection with pattern and calorimeter cylinder (r0=4000,
908  // c0=6000) if there is one, then track will be split
909 
910  bool intersect = MuonHoughMathUtils::lineThroughCylinder(patternpos, phi, theta, MuonHough::radius_cylinder, MuonHough::z_cylinder);
911 
912  if (!intersect) { // no split
913  ATH_MSG_DEBUG("Pattern not through calorimeter -> do not split ");
914  return {};
915  }
916 
917  return splitPatterns3D(phipattern, etapattern);
918 }
919 
920 std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::cleanupCombinedPattern(const Muon::MuonPrdPattern& combinedpattern) const {
921  const double phipattern = combinedpattern.globalDirection().phi();
922  const double thetapattern = combinedpattern.globalDirection().theta();
923 
924  CxxUtils::sincos scthetapattern(thetapattern);
925 
926  const Amg::Vector3D& patternpos = combinedpattern.globalPosition();
927  const double posx = patternpos.x();
928  const double posy = patternpos.y();
929  const double posz = patternpos.z();
930 
931  double invcurvature = 0.;
932  double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(posx, posy, phipattern);
933  double charge = 1.;
934  if (r0 < 0) charge = -1.;
935  double curvature = combinedpattern.globalDirection().mag();
936  if (curvature > 2) invcurvature = charge / curvature;
937 
938  ATH_MSG_DEBUG("cleaned up pattern: phi " << phipattern << " theta: " << thetapattern << " position: " << posx << " " << posy << " "
939  << posz);
940  ATH_MSG_DEBUG("Cleanup pattern charge " << charge << " curvature " << curvature);
941 
942  std::unique_ptr<Muon::MuonPrdPattern> combinedpattern_cleaned =
943  std::make_unique<Muon::MuonPrdPattern>(combinedpattern.globalPosition(), combinedpattern.globalDirection());
944 
945  for (unsigned int hitid = 0; hitid < combinedpattern.numberOfContainedPrds(); hitid++) {
946  const Trk::PrepRawData* prd = combinedpattern.prd(hitid);
947  const Amg::Vector3D& globalposhit = globalPos(prd);
948 
949  double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(posx, posy, phipattern);
950  double distance_xy = MuonHoughMathUtils::distanceToLine(globalposhit.x(), globalposhit.y(), r0, phipattern);
951 
952  double radius_pattern = globalposhit.perp();
953  double z0 = posz - radius_pattern * scthetapattern.cs / scthetapattern.sn;
954 
955  double distance_rz = MuonHoughMathUtils::signedDistanceCurvedToHit(z0, thetapattern, invcurvature, globalposhit);
956 
957  const double scale = std::max(1., globalposhit.mag() / 7000.);
958  ATH_MSG_VERBOSE("hit: " << globalposhit);
959  ATH_MSG_VERBOSE("CLEAN distancetopattern: "
960  << " dist xy " << distance_xy << " dist rz " << distance_rz << " scale: " << scale);
961  if (std::abs(distance_xy) < scale * m_maximum_xydistance && std::abs(distance_rz) < m_maximum_rzdistance) {
962  combinedpattern_cleaned->addPrd(prd);
963  } else {
964  ATH_MSG_DEBUG("Hit discarded: " << hitid << " dis xy " << distance_xy << " dis rz " << distance_rz);
965  ATH_MSG_DEBUG("Hit info: "<<m_idHelperSvc->toString(prd->identify()));
966  }
967  }
968 
969  ATH_MSG_DEBUG("size of cleaned pattern: " << combinedpattern_cleaned->numberOfContainedPrds());
970 
971  if (combinedpattern_cleaned->numberOfContainedPrds() == 0 && combinedpattern.numberOfContainedPrds() != 0) {
972  ATH_MSG_DEBUG(
973  "cleaned up pattern is empty (should happen only when initially no phi "
974  "pattern found and phi hits are added by ascociation map)");
975  }
976 
977  return combinedpattern_cleaned;
978 }
979 
980 std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::makeAssPhiPattern(
981  const Muon::MuonPrdPattern& muonpattern,
983  const EtaPhiHitAssocMap& phiEtaHitAssMap,
984  bool check_already_on_pattern) const {
985  // bool hits_added = false;
986  const unsigned int size = muonpattern.numberOfContainedPrds();
987 
988  PrepDataSet hits;
989  if (check_already_on_pattern) {
990  for (unsigned int i = 0; i < size; i++) { hits.insert(muonpattern.prd(i)); }
991  }
992  std::vector<const Trk::PrepRawData*> phihits;
993  for (unsigned int i = 0; i < size; i++) {
994  const Trk::PrepRawData* prd = muonpattern.prd(i);
995  // check on type of prd?
996  const Muon::MuonCluster* muoncluster = dynamic_cast<const Muon::MuonCluster*>(prd);
997  if (!muoncluster) continue;
998  EtaPhiHitAssocMap::const_iterator itr = phiEtaHitAssMap.find(prd);
1000  if (itr == phiEtaHitAssMap.end()) {
1001  continue;
1002  }
1003  std::copy_if(itr->second.begin(), itr->second.end(), std::back_inserter(phihits),
1004  [&check_already_on_pattern, &hits](const Trk::PrepRawData* phiHit){
1005  return !check_already_on_pattern || hits.insert(phiHit).second;
1006  });
1007  }
1008 
1009 
1010  if (phihits.empty()) { return nullptr; }
1011 
1012  std::unique_ptr<Muon::MuonPrdPattern> phipattern;
1013 
1014  double phi = 0., sin_phi = 0., cos_phi = 0.;
1015  if (m_use_cosmics) {
1016  phi = muonpattern.globalDirection().phi();
1017  } else {
1018  for (const Trk::PrepRawData* phiHit : phihits) {
1019  const Amg::Vector3D& globalposhit = globalPos(phiHit);
1020  CxxUtils::sincos scphihit(globalposhit.phi());
1021  sin_phi += scphihit.sn;
1022  cos_phi += scphihit.cs;
1023  }
1024  phi = std::atan2(sin_phi, cos_phi);
1025  }
1026 
1027  const double curvature = muonpattern.globalDirection().mag();
1028  const double theta = muonpattern.globalDirection().theta();
1029  CxxUtils::sincos scphi(phi);
1030  CxxUtils::sincos sctheta(theta);
1031 
1032  if (m_use_cosmics) {
1033  phipattern = std::make_unique<Muon::MuonPrdPattern>(muonpattern.globalPosition(), muonpattern.globalDirection(), phihits);
1034  } else {
1035  static const Amg::Vector3D globalpos{0.001, 0.001, 0.};
1036  const Amg::Vector3D& globaldir{curvature * scphi.cs * sctheta.sn,
1037  curvature * scphi.sn * sctheta.sn, curvature * sctheta.cs};
1038  phipattern = std::make_unique<Muon::MuonPrdPattern>(globalpos, globaldir, phihits);
1039  }
1040 
1041 
1042  // perform cleaning on newly created phipattern:
1043  std::unique_ptr<Muon::MuonPrdPattern> phipatternclean = cleanPhiPattern(std::move(phipattern));
1044 
1045  if (phipatternclean->numberOfContainedPrds() <= 0) {
1046  phipatternclean.reset();
1047  }
1048  return phipatternclean;
1049 }
1050 
1051 std::array<double,4> MuonCombinePatternTool::updateParametersForCosmics(const Muon::MuonPrdPattern& phipattern,
1052  const Muon::MuonPrdPattern& etapattern) const {
1053  // method returns r0, phi, rz0, theta
1054 
1055  const unsigned int etasize = etapattern.numberOfContainedPrds();
1056  const unsigned int phisize = phipattern.numberOfContainedPrds();
1057 
1058  const Amg::Vector3D& etaglobalpos = etapattern.globalPosition();
1059  const Amg::Vector3D& etaglobaldir = etapattern.globalDirection();
1060 
1061  const Amg::Vector3D& phiglobalpos = phipattern.globalPosition();
1062  const Amg::Vector3D& phiglobaldir = phipattern.globalDirection();
1063 
1064  std::array<double,4> old_pars{0};
1065 
1066  old_pars[0] = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phiglobaldir.phi());
1067  old_pars[1] = phiglobaldir.phi();
1068 
1069  const double theta_orig = etaglobaldir.theta();
1070  old_pars[2] = etaglobalpos.z() * std::sin(theta_orig);
1071  old_pars[3] = theta_orig;
1072 
1073  if (phisize + etasize <= 1) return old_pars;
1074 
1075  // first calculate new phi, r0 (initial value , -Pi/2):
1076  std::pair<double, double> rphi_start = calculateR0Phi(phipattern, etapattern);
1077  // for stabilising (can be cpu-optimised greatly):
1078  std::pair<double, double> rphi = calculateR0Phi(phipattern, etapattern, rphi_start.first);
1079 
1080  if (msgLvl(MSG::DEBUG) && std::abs(std::sin(rphi.first - rphi_start.first)) > 0.15 &&
1081  std::abs(std::sin(etaglobaldir.phi() - phiglobaldir.phi())) < 0.15) {
1082  ATH_MSG_DEBUG("unexpected phi change!");
1083  ATH_MSG_DEBUG("phi first: " << rphi_start.first << " phi second: " << rphi.first);
1084  ATH_MSG_DEBUG("phi etapattern: " << etaglobaldir.phi() << " phi phipattern: " << phiglobaldir.phi());
1085  }
1086 
1087  const double phi = rphi.first;
1088  const double r0 = rphi.second;
1089 
1090  CxxUtils::sincos scphi(phi);
1091 
1092  // calculate new theta and rz0: (not using phi hits this time, as eta pattern
1093  // is larger in general
1094 
1095  double av_radii{0.}, av_z{0.};
1096 
1097  for (unsigned int i = 0; i < etasize; ++i) {
1098  const Trk::PrepRawData* prd = etapattern.prd(i);
1099  const Amg::Vector3D globalposhit = globalPos(prd);
1100  av_radii += scphi.apply(globalposhit.y(), globalposhit.x());
1101  av_z += globalposhit.z();
1102  }
1103 
1104  if (etasize > 0) {
1105  av_radii /= etasize;
1106  av_z /= etasize;
1107  }
1108  double sumr = 0.;
1109  double sumz = 0.;
1110  for (unsigned int i = 0; i < etasize; i++) {
1111  const Trk::PrepRawData* prd = etapattern.prd(i);
1112  const Amg::Vector3D& globalposhit = globalPos(prd);
1113 
1114  double radius = scphi.apply(globalposhit.y(), globalposhit.x()); // hitx*scphi.cs + hity*scphi.sn;
1115  double r_offset = radius - av_radii;
1116  double z_offset = globalposhit.z() - av_z;
1117  double weight = r_offset * r_offset + z_offset * z_offset;
1118  int sign = 1;
1119  if (r_offset * std::cos(theta_orig) + z_offset * std::sin(theta_orig) < 0) { sign = -1; }
1120  sumr += weight * sign * r_offset;
1121  sumz += weight * sign * z_offset;
1122  }
1123 
1124  // const double sum_tan = sum_tanweight/sum_weight;
1125 
1126  ATH_MSG_VERBOSE("av_z : " << av_z << " av_radii: " << av_radii << " sumr: " << sumr << " sumz: " << sumz);
1127  if (std::abs(sumr) < 0.000001 || std::abs(sumz) < 0.000001) return old_pars;
1128 
1129  double theta = std::atan2(sumr, sumz);
1130 
1131  if (theta < 0) theta += M_PI;
1132 
1133  double rz0 = calculateRz0(etapattern, phi, theta);
1134 
1135  // ATH_MSG_DEBUG("old method rz0: " << sctheta.apply(av_z,-av_radii) );
1136  // const double rz0 = sctheta.apply(av_z,-av_radii); // (av_z * sctheta.sn) -
1137  // av_radii * sctheta.cs;
1138 
1139  std::array<double,4> new_pars {r0,phi,rz0,theta};
1140 
1141  ATH_MSG_VERBOSE("updated parameters: r0: " << new_pars[0] << " phi: " << new_pars[1] << " rz0: " << new_pars[2]
1142  << " theta: " << new_pars[3]);
1143  ATH_MSG_VERBOSE("old parameters: r0: " << old_pars[0] << " phi: " << old_pars[1] << " rz0: " << old_pars[2] << " theta: " << old_pars[3]);
1144 
1145  if (msgLvl(MSG::VERBOSE)) {
1146  ATH_MSG_VERBOSE("phisize: " << phisize << " etasize: " << etasize);
1147  for (unsigned int i = 0; i < phisize; i++) {
1148  const Trk::PrepRawData* prd = phipattern.prd(i);
1149  const Amg::Vector3D& globalposhit = globalPos(prd);
1150  double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), r0, phi);
1151  ATH_MSG_VERBOSE("distance to updated parameters in xy: " << distance);
1152  distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), old_pars[0], old_pars[1]);
1153  ATH_MSG_VERBOSE("old distance phi hit: " << distance);
1154  }
1155  for (unsigned int i = 0; i < etasize; i++) {
1156  const Trk::PrepRawData* prd = etapattern.prd(i);
1157  const Amg::Vector3D& globalposhit = globalPos(prd);
1158  double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), r0, phi);
1159  ATH_MSG_VERBOSE("distance to updated parameters in xy: " << distance);
1160  distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit.x(), globalposhit.y(), old_pars[0], old_pars[1]);
1161  ATH_MSG_VERBOSE("old distance eta hit: " << distance);
1162  }
1163  for (unsigned int i = 0; i < etasize; ++i) {
1164  const Trk::PrepRawData* prd = etapattern.prd(i);
1165  const Amg::Vector3D& globalposhit = globalPos(prd);
1166  double perp = scphi.apply(globalposhit.y(),
1167  globalposhit.x()); // globalposhit.x()*scphi.cs + globalposhit.y()*scphi.sn;
1168  double distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit[Amg::z], perp, rz0, theta);
1169  ATH_MSG_VERBOSE("distance to updated parameters in Rz: " << distance);
1170  distance = MuonHoughMathUtils::signedDistanceToLine(globalposhit[Amg::z], perp, old_pars[2], old_pars[3]);
1171  ATH_MSG_VERBOSE("old distance: " << distance);
1172  }
1173  }
1174  return new_pars;
1175 }
1176 
1177 std::pair<double, double> MuonCombinePatternTool::calculateR0Phi(const Muon::MuonPrdPattern& phipattern,
1178  const Muon::MuonPrdPattern& etapattern, double phi_est) const {
1179  // use eta pattern as well, since phi patterns consist sometimes of only 1
1180  // station etahit error 200mrad (2Pi/16*2), phi hit 20mrad (should be hough
1181  // binsize (18 mrad), but prefer ~ factor 10 for stabilility)
1182 
1183  // test if lever_arm > 2 m before updating , if not old values are used
1184 
1185  ATH_MSG_VERBOSE("calculateR0Phi");
1186 
1187  CxxUtils::sincos scphi_est(phi_est);
1188 
1189  const unsigned int etasize = etapattern.numberOfContainedPrds();
1190  const unsigned int phisize = phipattern.numberOfContainedPrds();
1191 
1192  const Amg::Vector3D& etaglobaldir = etapattern.globalDirection();
1193  const double phi_etapattern = etaglobaldir.phi();
1194  CxxUtils::sincos scphi_eta(phi_etapattern);
1195 
1196  const Amg::Vector3D& phiglobaldir = phipattern.globalDirection();
1197  const Amg::Vector3D& phiglobalpos = phipattern.globalPosition();
1198  const double phi_phipattern = phiglobaldir.phi();
1199  CxxUtils::sincos scphi_phi(phi_phipattern);
1200 
1201  const double phi_error_inv = 1. / 20.;
1202  const double phi_error_inv2 = phi_error_inv * phi_error_inv;
1203  const double eta_error_inv = 1. / 400.;
1204  const double eta_error_inv2 = eta_error_inv * eta_error_inv;
1205 
1206  // from MuonHoughPattern::updateParametersRPhi (partial code duplication.. :(
1207  // )
1208 
1209  double sum_etax{0.}, sum_etay{0.}, sum_phix{0.}, sum_phiy{0.};
1210 
1211  // calculate average point
1212 
1213  for (unsigned int i = 0; i < etasize; i++) {
1214  const Trk::PrepRawData* prd = etapattern.prd(i);
1215  const Amg::Vector3D& globalposhit = globalPos(prd);
1216  sum_etax += globalposhit.x();
1217  sum_etay += globalposhit.y();
1218  }
1219 
1220  for (unsigned int i = 0; i < phisize; i++) {
1221  const Trk::PrepRawData* prd = phipattern.prd(i);
1222  const Amg::Vector3D& globalposhit = globalPos(prd);
1223  sum_phix += globalposhit.x();
1224  sum_phiy += globalposhit.y();
1225  }
1226 
1227  const double av_x = (eta_error_inv2 * sum_etax + phi_error_inv2 * sum_phix) / (eta_error_inv2 * etasize + phi_error_inv2 * phisize);
1228  const double av_y = (eta_error_inv2 * sum_etay + phi_error_inv2 * sum_phiy) / (eta_error_inv2 * etasize + phi_error_inv2 * phisize);
1229 
1230  ATH_MSG_VERBOSE(" av_x: " << av_x << " av_y: " << av_y);
1231 
1232  // calculate weighted sum:
1233 
1234  double sumx {0.}, sumy {0.};
1235 
1236  // keep track of extreme points
1237 
1238  double x_min {0.}, x_max {0.}, y_min {0.}, y_max {0.}, lever_min {0.}, lever_max{0.};
1239 
1240  for (unsigned int i = 0; i < etasize; i++) {
1241  const Trk::PrepRawData* prd = etapattern.prd(i);
1242  const Amg::Vector3D& globalposhit = globalPos(prd);
1243  double x_offset = globalposhit.x() - av_x;
1244  double y_offset = globalposhit.y() - av_y;
1245  double height_squared = x_offset * x_offset + y_offset * y_offset;
1246  double weight = height_squared * eta_error_inv2;
1247  int sign = 1;
1248  if (x_offset * scphi_est.cs + y_offset * scphi_est.sn < 0) { sign = -1; }
1249  sumx += weight * sign * x_offset;
1250  sumy += weight * sign * y_offset;
1251 
1252  if (sign == 1 && height_squared > lever_max) {
1253  lever_max = height_squared;
1254  x_max = globalposhit.x();
1255  y_max = globalposhit.y();
1256  } else if (sign == -1 && height_squared > lever_min) {
1257  lever_min = height_squared;
1258  x_min = globalposhit.x();
1259  y_min = globalposhit.y();
1260  }
1261  }
1262 
1263  for (unsigned int i = 0; i < phisize; i++) {
1264  const Trk::PrepRawData* prd = phipattern.prd(i);
1265  const Amg::Vector3D& globalposhit = globalPos(prd);
1266  double x_offset = globalposhit.x() - av_x;
1267  double y_offset = globalposhit.y() - av_y;
1268  double height_squared = x_offset * x_offset + y_offset * y_offset;
1269  double weight = height_squared * phi_error_inv2;
1270  int sign = 1;
1271  if (x_offset * scphi_est.cs + y_offset * scphi_est.sn < 0) { sign = -1; }
1272  sumx += weight * sign * x_offset;
1273  sumy += weight * sign * y_offset;
1274 
1275  if (sign == 1 && height_squared > lever_max) {
1276  lever_max = height_squared;
1277  x_max = globalposhit.x();
1278  y_max = globalposhit.y();
1279  } else if (sign == -1 && height_squared > lever_min) {
1280  lever_min = height_squared;
1281  x_min = globalposhit.x();
1282  y_min = globalposhit.y();
1283  }
1284  }
1285 
1286  ATH_MSG_VERBOSE("av_x : " << av_x << " av_y: " << av_y << " sumx: " << sumx << " sumy: " << sumy);
1287 
1288  if (std::abs(sumx) < 0.000001 || std::abs(sumy) < 0.000001) {
1289  ATH_MSG_DEBUG(" sum too small to update");
1290 
1291  return std::make_pair(phi_phipattern, MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phi_phipattern));
1292  }
1293 
1294  // lever arm has to be larger than 2 m, else no update:
1295  if (std::hypot(x_max - x_min , y_max - y_min) < 2000) {
1296  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))
1297  << " x_max: " << x_max << " x_min: " << x_min << " y_max: " << y_max
1298  << " y_min: " << y_min);
1299  return std::make_pair(phi_phipattern, MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(phiglobalpos.x(), phiglobalpos.y(), phi_phipattern));
1300  }
1301 
1302  double phi_fit = std::atan2(sumy, sumx);
1303  if (phi_fit > 0) phi_fit -= M_PI; // phi between 0,-Pi for cosmics!
1304  CxxUtils::sincos scphi(phi_fit);
1305  const double r0_fit = scphi.apply(av_x, -av_y); // av_x * scphi.sn - av_y * scphi.cs;
1306 
1307  return std::make_pair(phi_fit, r0_fit);
1308 }
1309 
1310 double MuonCombinePatternTool::calculateRz0(const Muon::MuonPrdPattern& pattern, double phi, double theta) {
1311  double nhits = pattern.numberOfContainedPrds();
1312  CxxUtils::sincos sctheta(theta);
1313  CxxUtils::sincos scphi(phi);
1314 
1315  /*
1316  x = r_0 sin(phi) + t*cos(phi)sin(theta)
1317  y = - r_0 cos(phi) + t*sin(phi)sin(theta)
1318  z = z_0 + t*cos(theta)
1319 
1320  2 methods (average point of average_z0):
1321 
1322  r^2 = x^2+y^2 = r_0^2 + t^2*sin^2(theta)
1323  (this projects the hit radially onto the line)
1324 
1325  Not so good as the radius of the hits can be smaller than r_0
1326 
1327  method based on:
1328 
1329  x*cos(phi) + y*sin(phi) = t*sin(theta)
1330  (projects the hit on the line in x-y and calculates the distance to r_0)
1331 
1332  works best
1333  */
1334 
1335  // method 3:
1336 
1337  double rz0 = 0.;
1338  for (unsigned int i = 0; i < nhits; i++) {
1339  const Trk::PrepRawData* prd = pattern.prd(i);
1340  const Amg::Vector3D& poshit = globalPos(prd);
1341  int sign = (poshit.x() * scphi.cs + poshit.y() * scphi.sn > 0) ? 1 : -1;
1342  rz0 += poshit.z() * sctheta.sn - sign * sctheta.cs * poshit.perp();
1343  }
1344 
1345  if (nhits > 0) rz0 /= nhits;
1346  return rz0;
1347 }
1348 
1349 MuonCombinePatternTool::PrdPatternPair MuonCombinePatternTool::updatePatternsForCosmics(
1350  const Muon::MuonPrdPattern& phipattern, const Muon::MuonPrdPattern& etapattern, const std::array<double,4>& new_pars) const {
1351  double phi = new_pars[1];
1352  double theta = new_pars[3];
1353 
1354  CxxUtils::sincos scphi(phi);
1355  CxxUtils::sincos sctheta(theta);
1356 
1357  double x0 = new_pars[0] * scphi.sn;
1358  double y0 = -new_pars[0] * scphi.cs;
1359  double z0_phi = new_pars[2];
1360  double z0_eta = new_pars[2];
1361  if (std::abs(sctheta.sn) > 1e-7) {
1362  z0_phi = (new_pars[2] + new_pars[0] * sctheta.cs) / sctheta.sn; // z0 belonging to (x0,y0)
1363  z0_eta = new_pars[2] / sctheta.sn; // z0 of rz0
1364  }
1365 
1366  const Amg::Vector3D globalDir{scphi.cs * sctheta.sn, scphi.sn * sctheta.sn, sctheta.cs};
1367  const Amg::Vector3D globalPosPhi{x0, y0, z0_phi};
1368  const Amg::Vector3D globalPosEta{x0, y0, z0_eta};
1369  ATH_MSG_VERBOSE("updatePatternsForCosmics() -- eta pattern "<<std::endl
1370  <<m_printer->print(*etapattern.prd(0))<<std::endl
1371  <<" -- phi pattern "<<std::endl
1372  <<m_printer->print(*phipattern.prd(0))<<std::endl
1373  <<"Theta: " <<theta<<" globalPosEta: "<<Amg::toString(globalPosEta)<<", globalDir: "<<Amg::toString(globalDir));
1374 
1375 
1376  std::unique_ptr<Muon::MuonPrdPattern> updatedphipattern = std::make_unique<Muon::MuonPrdPattern>(globalPosPhi, globalDir, phipattern.prepRawDataVec());
1377  std::unique_ptr<Muon::MuonPrdPattern> updatedetapattern = std::make_unique<Muon::MuonPrdPattern>(globalPosEta, globalDir, etapattern.prepRawDataVec());
1378 
1379  return std::make_pair(std::move(updatedphipattern), std::move(updatedetapattern));
1380 }
1381 
1382 std::unique_ptr<MuonPatternCombinationCollection> MuonCombinePatternTool::makePatternCombinations(const MuonPrdPatternCollection& muonpatterns) const {
1383  ATH_MSG_DEBUG("makePatternCombinations");
1384 
1385  std::unique_ptr<MuonPatternCombinationCollection> patterncombinations = std::make_unique<MuonPatternCombinationCollection>();
1386 
1387  for (const Muon::MuonPattern* pit : muonpatterns) {
1388  const Amg::Vector3D roadmom = pit->globalDirection();
1389  const Amg::Vector3D roadpos = pit->globalPosition();
1390  ATH_MSG_DEBUG("phi: " << roadmom.phi() << " eta: " << roadmom.eta());
1391  ATH_MSG_DEBUG("x: " << roadpos.x() << " y: " << roadpos.y() << " z: " << roadpos.z());
1392 
1393  // sort pattern per chamber
1394  std::map<Identifier, std::vector<const Trk::PrepRawData*>> chamberMap;
1395  for (unsigned int i = 0; i < pit->numberOfContainedPrds(); ++i) {
1396  const Trk::PrepRawData* prd = pit->prd(i);
1397  Identifier channelId = prd->identify();
1398  const Identifier moduleId = m_idHelperSvc->chamberId(channelId);
1399  std::vector<const Trk::PrepRawData*>& chambVec = chamberMap[moduleId];
1400  if (chambVec.empty()) chambVec.reserve(pit->numberOfContainedPrds());
1401  chambVec.push_back(prd);
1402  }
1403 
1404  // build chamberintersect vector
1405  std::vector<Muon::MuonPatternChamberIntersect> mpciVec;
1406  mpciVec.reserve(chamberMap.size());
1407  Amg::Vector3D patpose{Amg::Vector3D::Zero()}, patdire{Amg::Vector3D::Zero()};
1408  for ( const auto& [moduleId, chambVec] : chamberMap) {
1409  const Trk::PrepRawData* prd = chambVec.front();
1410  const Amg::Vector3D& globalpos = globalPos(prd);
1411  if (m_use_cosmics) {
1412  // not flip
1413  patdire = roadmom.unit();
1414  patpose = roadpos;
1415  } else {
1416  MuonHoughMathUtils::extrapolateCurvedRoad(roadpos, roadmom, globalpos, patpose, patdire);
1417  }
1418 
1419  Muon::MuonPatternChamberIntersect mpci = Muon::MuonPatternChamberIntersect(patpose, patdire, chambVec);
1420  mpciVec.push_back(mpci);
1421  }
1422 
1423  ATH_MSG_VERBOSE(__FILE__<<":"<<__LINE__<<" Pattern position "<<Amg::toString(roadpos));
1424  Amg::Vector3D tvertex{Amg::Vector3D::Zero()};
1425  Trk::TrackParameters* parameters = new Trk::Perigee(roadpos, roadmom, 1., tvertex); // if -1 then charge flipped anyway
1426  Muon::MuonPatternCombination* combination = new Muon::MuonPatternCombination(parameters, mpciVec);
1427  patterncombinations->push_back(combination);
1428  }
1429  return patterncombinations;
1430 }
1431 
1432 bool MuonCombinePatternTool::subset(const Muon::MuonPrdPattern* pattern1, const Muon::MuonPrdPattern* pattern2) {
1435  // first check if pattern 1 is not larger than 2:
1436  if (pattern1->numberOfContainedPrds() > pattern2->numberOfContainedPrds()) { return false; }
1437 
1438  PrepDataSet hits;
1439  for (unsigned int hitnr = 0; hitnr < pattern2->numberOfContainedPrds(); ++hitnr) { hits.insert(pattern2->prd(hitnr)); }
1440 
1441  for (unsigned int hitnr = 0; hitnr < pattern1->numberOfContainedPrds(); ++hitnr) {
1442  if (!hits.count(pattern1->prd(hitnr))) {
1443  return false;
1444  }
1445  }
1446  return true;
1447 }
1448 
1449 bool MuonCombinePatternTool::subset(std::pair<PrepDataSet,
1450  PrepDataSet>& candidate1,
1451  std::pair<PrepDataSet,
1452  PrepDataSet>& candidate2) {
1455  // first check if pattern 1 is not larger than 2:
1456  if (candidate1.first.size() > candidate2.first.size() ||
1457  candidate1.second.size() > candidate2.second.size()) { return false; }
1458  for (const Trk::PrepRawData* find_me : candidate1.first) {
1459  if (!candidate2.first.count(find_me)) { return false; }
1460  }
1461  for (const Trk::PrepRawData* find_me : candidate1.second) {
1462  if (candidate2.second.count(find_me)) { return false; }
1463  }
1464 
1465  return true;
1466 }
1467 
1468 void MuonCombinePatternTool::printPattern(const Muon::MuonPrdPattern* muonpattern) const {
1469  if (msgLvl(MSG::VERBOSE)) {
1470  ATH_MSG_VERBOSE("Printout of Pattern: ");
1471 
1472  const Amg::Vector3D& pos = muonpattern->globalPosition();
1473  const Amg::Vector3D& dir = muonpattern->globalDirection();
1474 
1475  ATH_MSG_VERBOSE("pos: x: " << pos.x() << " y: " << pos.y() << " z: " << pos.z());
1476  ATH_MSG_VERBOSE("dir: x: " << dir.x() << " y: " << dir.y() << " z: " << dir.z());
1477  ATH_MSG_VERBOSE("phi: " << dir.phi() << " theta: " << dir.theta() << " rz0: " << pos.z() * std::sin(dir.theta()));
1478 
1479  for (unsigned int k = 0; k < muonpattern->numberOfContainedPrds(); k++) {
1480  const Trk::PrepRawData* prd = muonpattern->prd(k);
1481  const Muon::MdtPrepData* mdtprd = dynamic_cast<const Muon::MdtPrepData*>(prd);
1482  if (mdtprd) {
1483  const Trk::Surface& surface = mdtprd->detectorElement()->surface(mdtprd->identify());
1484  const Amg::Vector3D& gpos = surface.center();
1485  ATH_MSG_VERBOSE("mdt " << k << " x: " << gpos.x() << " y: " << gpos.y() << " z: " << gpos.z());
1486  } else if (!mdtprd) {
1487  const Muon::MuonCluster* muoncluster = dynamic_cast<const Muon::MuonCluster*>(prd);
1488  if (muoncluster) {
1489  const Amg::Vector3D& gpos = muoncluster->globalPosition();
1490  ATH_MSG_VERBOSE("cluster " << k << " x: " << gpos.x() << " y: " << gpos.y() << " z: " << gpos.z());
1491  }
1492  if (!muoncluster) { ATH_MSG_VERBOSE("no muon prd? "); }
1493  }
1494  }
1495  }
1496 }
1497 
1498 std::unique_ptr<Muon::MuonPrdPattern> MuonCombinePatternTool::cleanPhiPattern(std::unique_ptr<Muon::MuonPrdPattern> phipattern) const {
1499  const Amg::Vector3D& olddir = phipattern->globalDirection();
1500  const double theta = olddir.theta();
1501  const unsigned int size = phipattern->numberOfContainedPrds();
1502 
1503  if (msgLvl(MSG::DEBUG)) {
1504  ATH_MSG_DEBUG("Start Phi hits cleaning with " << size << " hits "
1505  << " theta " << theta);
1506  const Amg::Vector3D& oldpos = phipattern->globalPosition();
1507  double r0 = MuonHoughMathUtils::signedDistanceOfLineToOrigin2D(oldpos.x(), oldpos.y(), olddir.phi());
1508  ATH_MSG_DEBUG("Start Phi: " << olddir.phi() << " r0: " << r0);
1509  }
1510 
1511  // need internal class to be able to remove hits fast
1512  std::unique_ptr<MuonHoughPattern> newpattern = std::make_unique<MuonHoughPattern>(MuonHough::hough_xy);
1513  for (unsigned int phihitnr = 0; phihitnr < size; phihitnr++) { newpattern->addHit(std::make_shared<MuonHoughHit>(phipattern->prd(phihitnr))); }
1514  newpattern->updateParametersRPhi(m_use_cosmics);
1515  double phi = newpattern->getEPhi();
1516  double r0 = newpattern->getERPhi();
1517 
1518  CxxUtils::sincos scphi(phi);
1519 
1520  constexpr int number_of_iterations = 4;
1521  std::array<double,number_of_iterations> cutvalues {1000., 500., 250., 125.};
1522  if (m_use_cosmics) {
1523  cutvalues[0] = 5000.;
1524  cutvalues[1] = 2500.;
1525  cutvalues[2] = 1250.;
1526  cutvalues[3] = 1250.;
1527  }
1528 
1529  for (int it = 0; it < number_of_iterations; ++it) {
1530  ATH_MSG_VERBOSE("iteration " << it << " cutvalue: " << cutvalues[it]);
1531  bool change = true;
1532  while (change) {
1533  ATH_MSG_VERBOSE("size: " << newpattern->size() << " r0: " << r0 << " phi: " << phi);
1534 
1535  double max_dist = 0.;
1536  unsigned int max_i = 99999;
1537  for (unsigned int i = 0; i < newpattern->size(); i++) {
1538  double dist =
1539  scphi.apply(newpattern->getHitx(i), -newpattern->getHity(i)) - r0; // newpattern->getHitx(i) * scphi.sn -
1540  // newpattern->getHity(i) * scphi.cs - r0;
1541  ATH_MSG_VERBOSE("Dist: " << dist);
1542  if (std::fabs(dist) > std::abs(max_dist)) {
1543  max_dist = dist;
1544  max_i = i;
1545  }
1546  }
1547  if (std::abs(max_dist) < cutvalues[it]) {
1548  change = false;
1549  } else {
1550  newpattern->removeHit(max_i);
1551  newpattern->updateParametersRPhi(m_use_cosmics);
1552  phi = newpattern->getEPhi();
1553  r0 = newpattern->getERPhi();
1554  scphi = CxxUtils::sincos(phi);
1555  }
1556  }
1557  }
1558 
1559  ATH_MSG_DEBUG("Final size: " << newpattern->size() << " r0: " << r0 << " phi: " << phi);
1560 
1561  // update parameters rz (not very important as later overwritten)
1562  // put r0 to IP for collisions
1563 
1564  double thetanew = 0.;
1565  double r0_new = 1.; // put 1 mm r0 value
1566 
1567  if (m_use_cosmics) { r0_new = r0; }
1568 
1569  unsigned int nPatterns = newpattern->size();
1570  for (unsigned int i = 0; i < nPatterns; i++) { thetanew += newpattern->getTheta(i); }
1571 
1572  if (nPatterns > 0) thetanew /= nPatterns;
1573 
1574  double z0_new = 0.;
1575 
1576  double x0_new = r0_new * scphi.sn;
1577  double y0_new = -r0_new * scphi.cs;
1578  CxxUtils::sincos sctheta(thetanew);
1579 
1580  const Amg::Vector3D pos {x0_new, y0_new, z0_new};
1581  const Amg::Vector3D dir {sctheta.sn * scphi.cs, sctheta.sn * scphi.sn, sctheta.cs};
1582 
1583  std::unique_ptr<Muon::MuonPrdPattern> cleanpattern = std::make_unique<Muon::MuonPrdPattern>(pos, dir);
1584 
1585  for (unsigned int i = 0; i < newpattern->size(); i++) { cleanpattern->addPrd(newpattern->getPrd(i)); }
1586 
1587  return cleanpattern;
1588 }
1589 
1590 void MuonCombinePatternTool::addCandidate(const CandPrdPatPtr& etapattern, const CandPrdPatPtr& phipattern,
1591  std::vector<CandidatePatPair>& candidates, bool add_asspattern,
1592  const EtaPhiHitAssocMap& phiEtaHitAssMap) const {
1593  if (!m_use_cosmics || !m_splitpatterns) {
1594  candidates.emplace_back(etapattern, phipattern);
1595  return;
1596  }
1597 
1598  std::vector<PrdPatternPair> splitpatterns = splitPatternsCylinder(phipattern.get(), etapattern.get());
1599 
1600  if (splitpatterns.empty()) {
1601  candidates.emplace_back(etapattern, phipattern);
1602  }
1603 
1604  else {
1605  for (auto& [phiPattern, etaPattern] : splitpatterns) {
1606  // skip when empty eta pattern , possible duplication when associated phi
1607  // pattern is found, but then will be cleaned later
1608  if (etaPattern->numberOfContainedPrds() == 0) {
1609  continue;
1610  }
1611  candidates.emplace_back(std::move(etaPattern), std::move(phiPattern));
1612  }
1613  }
1614 
1615  // make associated pattern don't split eta pattern yet, but split based on phi
1616  // of ass. pattern bool asspattern_added = false;
1617  if (!add_asspattern) { return ;}
1618  std::unique_ptr<Muon::MuonPrdPattern> assphipattern = makeAssPhiPattern(*etapattern, phiEtaHitAssMap, true);
1619 
1620  if (!assphipattern) { return;}
1621 
1622  // print associated pattern:
1623  if (msgLvl(MSG::VERBOSE)) {
1624  ATH_MSG_VERBOSE("Associated Pattern: ");
1625  printPattern(assphipattern.get());
1626  }
1627 
1628  // update parameters:
1629  std::array<double,4> new_pars = updateParametersForCosmics(*assphipattern, *etapattern);
1630  PrdPatternPair updatedpatterns = updatePatternsForCosmics(*assphipattern, *etapattern, new_pars);
1631 
1632  std::unique_ptr<Muon::MuonPrdPattern>& cosmicPhiPattern = updatedpatterns.first;
1633  std::unique_ptr<Muon::MuonPrdPattern>& cosmicEtaPattern = updatedpatterns.second;
1634 
1635  std::vector<PrdPatternPair> splitpatterns_ass = splitPatternsCylinder(cosmicPhiPattern.get(), cosmicEtaPattern.get());
1636  if (splitpatterns_ass.empty()) {
1637  candidates.emplace_back(std::move(cosmicEtaPattern), std::move(cosmicPhiPattern));
1638  return;
1639  }
1640 
1641 
1642  for (auto& [splitPhiPattern, splitEtaPattern] : splitpatterns_ass) {
1643  if (splitPhiPattern->numberOfContainedPrds() == 0 ||
1644  splitEtaPattern->numberOfContainedPrds() == 0) {
1645  continue;
1646  }
1647  candidates.emplace_back(std::move(splitEtaPattern), std::move(splitPhiPattern));
1648  }
1649 }
1650 
1651 void MuonCombinePatternTool::cleanCandidates(std::vector<CandidatePatPair>& candidates) {
1652  std::vector<CandidatePatPair>::iterator it1;
1653  std::vector<CandidatePatPair>::iterator it2;
1654  // map between set of prd's (eta and phi) and candidates , stored for speed
1655  std::map<CandidatePatPair, std::pair<PrepDataSet, PrepDataSet>> hitsMap;
1656 
1657  // fill map
1658  for (it1 = candidates.begin(); it1 != candidates.end(); ++it1) {
1659  PrepDataSet etahits;
1660  for (unsigned int hitnr = 0; hitnr < (*it1).first->numberOfContainedPrds(); hitnr++) { etahits.insert((*it1).first->prd(hitnr)); }
1661  PrepDataSet phihits;
1662  if ((*it1).second) { // phi pattern might be 0!
1663  for (unsigned int hitnr = 0; hitnr < (*it1).second->numberOfContainedPrds(); hitnr++) {
1664  phihits.insert((*it1).second->prd(hitnr));
1665  }
1666  }
1667  hitsMap.insert(std::make_pair((*it1), std::make_pair(etahits, phihits)));
1668  }
1669  // cppcheck-suppress invalidContainer; false positive
1670  for (it1 = candidates.begin(); it1 != candidates.end(); ++it1) {
1671  std::pair<PrepDataSet, PrepDataSet>& hits1 = hitsMap[(*it1)];
1672  it2 = it1 + 1;
1673  while (it2 != candidates.end()) {
1674  std::pair<PrepDataSet, PrepDataSet>& hits2 = hitsMap[(*it2)];
1675  if (subset((hits2), (hits1))) { // 2 subset of 1, remove 2 // in case of
1676  // equality best (earliest) is kept
1677  it2 = candidates.erase(it2); // it2 points to next item, it1 not invalidated!
1678  } else if (subset((hits1), (hits2))) { // 1 subset of 2, remove 1
1679  it1 = candidates.erase(it1); // it1 points to next item, it2 invalidated!
1680  it2 = it1 + 1;
1681  //cppcheck-suppress selfAssignment
1682  hits1 = hitsMap[(*it1)]; // redefine hits1
1683  } else {
1684  ++it2;
1685  }
1686  }
1687  }
1688 }
CxxUtils::sincos::apply
double apply(double a, double b) const
Definition: sincos.h:98
MuonCombinePatternTool::ChamberInfo::neta
int neta
Definition: MuonCombinePatternTool.h:21
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
mergePhysValFiles.pattern
pattern
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:26
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
Muon::MuonPattern::globalDirection
const Amg::Vector3D & globalDirection() const
Global direction of the pattern.
Definition: MuonPattern.h:58
Trk::PrepRawDataType::MdtPrepData
@ MdtPrepData
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
MuonCombinePatternTool::ChamberInfo::phiMax
double phiMax
Definition: MuonCombinePatternTool.h:28
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
MuonCombinePatternTool.h
MuonCombinePatternTool::m_maxSizePhiPatternLoose
unsigned int m_maxSizePhiPatternLoose
Definition: MuonCombinePatternTool.h:158
MuonCombinePatternTool::CandPrdPatPtr
std::shared_ptr< const Muon::MuonPrdPattern > CandPrdPatPtr
Definition: MuonCombinePatternTool.h:35
MuonPatternChamberIntersect.h
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
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
MuonHoughHitContainer::getTheta
double getTheta(unsigned int hitno) const
returns theta of hit hitno
Definition: MuonHoughHitContainer.h:99
max
#define max(a, b)
Definition: cfImp.cxx:41
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
TrackParameters.h
MuonCombinePatternTool::m_maximum_xydistance
const double m_maximum_xydistance
distance cut in xy for hits
Definition: MuonCombinePatternTool.h:135
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
Surface.h
Muon::MuonPrdPattern::addPrd
virtual void addPrd(const Trk::PrepRawData *prd)
add hit to pattern
Definition: MuonPrdPattern.h:60
perp
Scalar perp() const
perp method - perpenticular length
Definition: AmgMatrixBasePlugin.h:44
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
drawFromPickle.candidates
candidates
Definition: drawFromPickle.py:271
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
MuonHoughPattern::getEPhi
double getEPhi() const
returns phi of pattern
Definition: MuonHoughPattern.h:116
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
skel.it
it
Definition: skel.GENtoEVGEN.py:396
M_PI
#define M_PI
Definition: ActiveFraction.h:11
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
MuonHough::z_cylinder
constexpr double z_cylinder
length of cylinder
Definition: MuonHoughMathUtils.h:26
Trk::PrepRawData::type
virtual bool type(PrepRawDataType type) const =0
Interface method checking the type.
MuonCombinePatternTool::m_maximum_rzdistance
const double m_maximum_rzdistance
distance cut in rz for hits
Definition: MuonCombinePatternTool.h:137
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
MuonCombinePatternTool::initialize
StatusCode initialize() override
Definition: MuonCombinePatternTool.cxx:58
Muon::MuonCluster::globalPosition
virtual const Amg::Vector3D & globalPosition() const =0
Returns the global position of the measurement (calculated on the fly)
MuonPrepDataContainer.h
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
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
PrdPatternPair
MuonCombinePatternTool::PrdPatternPair PrdPatternPair
Definition: MuonCombinePatternTool.cxx:36
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:1051
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:1590
Muon::MuonPrdPattern
Class to store a pattern in the muon system containing PrepRawData pointers.
Definition: MuonPrdPattern.h:27
MuonCombinePatternTool::cleanCandidates
static void cleanCandidates(std::vector< CandidatePatPair > &candidates)
clean candidates from subsets or duplicates
Definition: MuonCombinePatternTool.cxx:1651
MuonHoughHitContainer::addHit
void addHit(const std::shared_ptr< MuonHoughHit > &hit)
add hit to container
Definition: MuonHoughHitContainer.cxx:8
MuonCombinePatternTool::subset
static bool subset(const Muon::MuonPrdPattern *pattern1, const Muon::MuonPrdPattern *pattern2)
is pattern1 a complete subset of other pattern2?
Definition: MuonCombinePatternTool.cxx:1432
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
MuonGM::MdtReadoutElement::getWireLength
double getWireLength(const int tubeLayer, const int tube) const
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
Muon::MuonPattern
Basic class for patterns in the muon spectrometer consistig out of a list of Trk::PrepRawData objects...
Definition: MuonPattern.h:23
Muon::IMuonCombinePatternTool::EtaPhiHitAssocMap
std::map< const Trk::PrepRawData *, PrepDataSet > EtaPhiHitAssocMap
Definition: IMuonCombinePatternTool.h:25
Amg::z
@ z
Definition: GeoPrimitives.h:36
MuonHough::hough_xy
@ hough_xy
Definition: MuonHoughPattern.h:14
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:1349
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
python.SystemOfUnits.meter
int meter
Definition: SystemOfUnits.py:61
CxxUtils::sincos::cs
double cs
Definition: sincos.h:95
EventPrimitivesToStringConverter.h
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
MuonGM::MdtReadoutElement
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MdtReadoutElement.h:51
MuonCombinePatternTool::m_use_cosmics
bool m_use_cosmics
use cosmic settings
Definition: MuonCombinePatternTool.h:140
lumiFormat.i
int i
Definition: lumiFormat.py:85
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
Muon::MuonPatternChamberIntersect
This class holds information needed for the Moore and MoMu pattern recognition for a muon chamber.
Definition: MuonPatternChamberIntersect.h:38
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
DetDescrDictionaryDict::it1
std::vector< HWIdentifier >::iterator it1
Definition: DetDescrDictionaryDict.h:17
MuonHoughHitContainer::size
unsigned int size() const
returns size of hitcontainer
Definition: MuonHoughHitContainer.h:104
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
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
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
MuonHoughHitContainer::getHity
double getHity(unsigned int hitno) const
returns y position of hit hitno
Definition: MuonHoughHitContainer.h:95
MuonGM::MdtReadoutElement::getActiveTubeLength
double getActiveTubeLength(const int tubeLayer, const int tube) const
MuonHoughHitContainer::getPrd
const Trk::PrepRawData * getPrd(unsigned int hitno) const
returns preprawdata pointer of hit hitno
Definition: MuonHoughHitContainer.h:110
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition: TRT_PAI_physicsConstants.h:20
test_pyathena.parent
parent
Definition: test_pyathena.py:15
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
MuonCombinePatternTool::PrdPatternPair
std::pair< std::unique_ptr< Muon::MuonPrdPattern >, std::unique_ptr< Muon::MuonPrdPattern > > PrdPatternPair
Definition: MuonCombinePatternTool.h:33
AnalysisUtils::copy_if
Out copy_if(In first, const In &last, Out res, const Pred &p)
Definition: IFilterUtils.h:30
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ParametersBase
Definition: ParametersBase.h:55
MuonCombinePatternTool::m_flipdirectionforcosmics
bool m_flipdirectionforcosmics
flip direction for cosmics after splitting as if coming from IP (false by default)
Definition: MuonCombinePatternTool.h:155
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
Muon::MuonPattern::globalPosition
const Amg::Vector3D & globalPosition() const
Global position of the pattern.
Definition: MuonPattern.h:57
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
CaloCondBlobAlgs_fillNoiseFromASCII.channelId
channelId
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:122
MuonHoughPattern::updateParametersRPhi
void updateParametersRPhi(bool cosmics=false)
update parameters in rphi plane based on weighted fit
Definition: MuonHoughPattern.cxx:192
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
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:778
beamspotman.dir
string dir
Definition: beamspotman.py:623
min
#define min(a, b)
Definition: cfImp.cxx:40
Trk::PrepRawData
Definition: PrepRawData.h:62
MuonCombinePatternTool::printPattern
void printPattern(const Muon::MuonPrdPattern *muonpattern) const
print out pattern hits
Definition: MuonCombinePatternTool.cxx:1468
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
MuonCombinePatternTool::cleanupCombinedPattern
std::unique_ptr< Muon::MuonPrdPattern > cleanupCombinedPattern(const Muon::MuonPrdPattern &combinedpattern) const
clean combined pattern, remove outliers
Definition: MuonCombinePatternTool.cxx:920
MuonCombinePatternTool::ChamberInfo::phiMin
double phiMin
Definition: MuonCombinePatternTool.h:27
MuonCombinePatternTool::ChamberInfo::ninsidePat
int ninsidePat
Definition: MuonCombinePatternTool.h:25
MuonCombinePatternTool::ChamberInfo::nphi
int nphi
Definition: MuonCombinePatternTool.h:22
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
MuonCombinePatternTool::m_splitpatterns
bool m_splitpatterns
split patterns (only for cosmics)
Definition: MuonCombinePatternTool.h:143
charge
double charge(const T &p)
Definition: AtlasPID.h:538
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:1498
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
Muon::MdtPrepData
Class to represent measurements from the Monitored Drift Tubes.
Definition: MdtPrepData.h:33
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MuonCombinePatternTool::MuonCombinePatternTool
MuonCombinePatternTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: MuonCombinePatternTool.cxx:38
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
MuonDetectorManager.h
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
MuonGM::MdtReadoutElement::ROPos
Amg::Vector3D ROPos(const int tubelayer, const int tube) const
Definition: MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:289
CxxUtils::sincos::sn
double sn
Definition: sincos.h:92
MuonCombinePatternTool::m_printer
PublicToolHandle< Muon::MuonEDMPrinterTool > m_printer
Definition: MuonCombinePatternTool.h:163
MuonCombinePatternTool::ChamberInfo::ninside
int ninside
Definition: MuonCombinePatternTool.h:23
Trk::inside
@ inside
Definition: PropDirection.h:29
MuonCombinePatternTool::m_useTightAssociation
bool m_useTightAssociation
Definition: MuonCombinePatternTool.h:157
Muon::MuonPrdPattern::PrdVector
std::vector< const Trk::PrepRawData * > PrdVector
Definition: MuonPrdPattern.h:29
MuonCombinePatternTool::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MuonCombinePatternTool.h:161
MuonCombinePatternTool::ChamberInfo::noutsidePat
int noutsidePat
Definition: MuonCombinePatternTool.h:26
MuonGM::MdtReadoutElement::tubePos
Amg::Vector3D tubePos(const Identifier &id) const
Returns the global position of the given tube.
MuonHoughPattern.h
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 of closest approach of two lines.
Definition: GeoPrimitivesHelpers.h:325
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
MuonCombinePatternTool::calculateRz0
static double calculateRz0(const Muon::MuonPrdPattern &pattern, double phi, double theta)
calculate rz0 for cosmic pattern
Definition: MuonCombinePatternTool.cxx:1310
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
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
MuonHoughHitContainer::removeHit
void removeHit(unsigned int hitno)
remove hit from container
Definition: MuonHoughHitContainer.cxx:13
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:1177
DEBUG
#define DEBUG
Definition: page_access.h:11
MuonCombinePatternTool::m_nodiscarding
bool m_nodiscarding
don't discard any candidates based on gasgap assocation (true by default)
Definition: MuonCombinePatternTool.h:147
MuonGM::MdtReadoutElement::surface
virtual const Trk::Surface & surface() const override final
Return surface associated with this detector element.
Definition: MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:891
MuonHoughPattern::getERPhi
double getERPhi() const
returns r0/d0 of pattern
Definition: MuonHoughPattern.h:117
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:76
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:889
Muon::MuonPrdPattern::numberOfContainedPrds
virtual unsigned int numberOfContainedPrds() const
Number or PrepRawData contained by this Pattern.
Definition: MuonPrdPattern.h:64
MuonCombinePatternTool::makePatternCombinations
virtual std::unique_ptr< MuonPatternCombinationCollection > makePatternCombinations(const MuonPrdPatternCollection &muonpatterns) const override
converts MuonPrdPatterns into MuonPatternCombinationCollection MuonPatternCombinationCollection are d...
Definition: MuonCombinePatternTool.cxx:1382
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Muon::IMuonCombinePatternTool::PrepDataSet
std::set< const Trk::PrepRawData *, IdentifierPrdLess > PrepDataSet
Definition: IMuonCombinePatternTool.h:24
Muon::MuonStationIndex::StIndex
StIndex
enum to classify the different station layers in the muon spectrometer
Definition: MuonStationIndex.h:23
Muon::MuonCluster
Class representing clusters in the muon system.
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPrepRawData/MuonPrepRawData/MuonCluster.h:37
MuonHough::radius_cylinder
constexpr double radius_cylinder
radius of cylinder
Definition: MuonHoughMathUtils.h:24
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
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
Muon::MuonPatternCombination
The MuonPatternCombination class provides the means to store the output of the initial global pattern...
Definition: MuonPatternCombination.h:29
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
MuonHoughHitContainer::getHitx
double getHitx(unsigned int hitno) const
returns x position of hit hitno
Definition: MuonHoughHitContainer.h:94
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
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:980
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
Muon::MdtPrepData::detectorElement
virtual const MuonGM::MdtReadoutElement * detectorElement() const override
Returns the detector element corresponding to this PRD.
Definition: MdtPrepData.h:141
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
calibdata.tube
tube
Definition: calibdata.py:31
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.
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
MuonCombinePatternTool::ChamberInfo::noutside
int noutside
Definition: MuonCombinePatternTool.h:24
Muon::MuonPrdPattern::prepRawDataVec
const PrdVector & prepRawDataVec() const
Definition: MuonPrdPattern.h:73
MuonCombinePatternTool::m_maxSizeEtaPatternLoose
unsigned int m_maxSizeEtaPatternLoose
Definition: MuonCombinePatternTool.h:159
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
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::ChamberInfo
Definition: MuonCombinePatternTool.h:19
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...
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Thu Nov 7 2024 21:21:52 for ATLAS Offline Software by   doxygen 1.8.18