ATLAS Offline Software
MSVertexTrackletTool.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 "MSVertexTrackletTool.h"
6 
7 #include "TMath.h"
8 #include "CxxUtils/trapping_fp.h"
9 
10 /*
11  Tracklet reconstruction tool
12  See documentation at https://cds.cern.ch/record/1455664 and https://cds.cern.ch/record/1520894
13 */
14 
15 namespace Muon {
16 
17  //** ----------------------------------------------------------------------------------------------------------------- **//
18 
19  // Delta Alpha Constants -- p = k/(delta_alpha)
20  // inner and outer small stations don't have a from data determined constant. Instead a default value for DeltaAlpgaCut, momentum and momentum error is used.
21  constexpr double c_BIL = 28.4366; // MeV*mrad
22  constexpr double c_BMS = 53.1259; // MeV*mrad
23  constexpr double c_BML = 62.8267; // MeV*mrad
24  constexpr double c_BOL = 29.7554; // MeV*mrad
25  //** ----------------------------------------------------------------------------------------------------------------- **//
26 
27  MSVertexTrackletTool::MSVertexTrackletTool(const std::string& type, const std::string& name, const IInterface* parent) :
28  AthAlgTool(type, name, parent) {
29  declareInterface<IMSVertexTrackletTool>(this);
30  }
31 
32  //** ----------------------------------------------------------------------------------------------------------------- **//
33 
34  StatusCode MSVertexTrackletTool::initialize() {
35  ATH_CHECK(m_mdtTESKey.initialize());
36  ATH_CHECK(m_TPContainer.initialize());
37  ATH_CHECK(m_idHelperSvc.retrieve());
38 
39  return StatusCode::SUCCESS;
40  }
41 
42  //** ----------------------------------------------------------------------------------------------------------------- **//
43 
44  StatusCode MSVertexTrackletTool::findTracklets(std::vector<Tracklet>& tracklets, const EventContext& ctx) const {
45  // record TrackParticle container in StoreGate
46  SG::WriteHandle<xAOD::TrackParticleContainer> container(m_TPContainer, ctx);
47  ATH_CHECK(container.record(std::make_unique<xAOD::TrackParticleContainer>(), std::make_unique<xAOD::TrackParticleAuxContainer>()));
48 
49  // sort the MDT hits into chambers & MLs
50  std::vector<std::vector<const Muon::MdtPrepData*> > SortedMdt;
51 
52  int nMDT = SortMDThits(SortedMdt, ctx);
53 
54  if (nMDT <= 0) { return StatusCode::SUCCESS; }
55 
56  if (msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG("MDT hits are selected and sorted");
57 
58  // loop over the MDT hits and find segments
59  // select the tube combinations to be fit
60  /*Select hits in at least 2 layers and require hits be ordered by increasing tube number (see diagrams below).
61  ( )( )(3)( ) ( )(3)( )( ) ( )( )( )( ) ( )( )(3)( ) ( )(2)(3)( )
62  ( )(2)( )( ) ( )(2)( )( ) ( )(2)(3)( ) ( )(1)(2)( ) ( )(1)( )( )
63  (1)( )( )( ) (1)( )( )( ) (1)( )( )( ) ( )( )( )( ) ( )( )( )( )
64  Barrel selection criteria: |z_mdt1 - z_mdt2| < m_d12_max (50 mm), |z_mdt1 - z_mdt3| < m_d13_max (80 mm)
65  Endcap selection criteria: |r_mdt1 - r_mdt2| < m_d12_max (50 mm), |r_mdt1 - r_mdt3| < m_d13_max (80 mm)
66  */
67 
68  std::vector<TrackletSegment> segs[6][2][16]; // single ML segment array (indicies [station type][ML][sector]) with station type iterating through barrel inner, middle, outer and then endcap inner, middle, outer
69  std::vector<std::vector<const Muon::MdtPrepData*> >::const_iterator ChamberItr = SortedMdt.begin();
70  for (; ChamberItr != SortedMdt.end(); ++ChamberItr) {
71  std::vector<TrackletSegment> mlsegments;
72  std::vector<const Muon::MdtPrepData*>::const_iterator mdt1 = ChamberItr->begin();
73  std::vector<const Muon::MdtPrepData*>::const_iterator mdtEnd = ChamberItr->end();
74  if (IgnoreMDTChamber(*mdt1)) continue;
75 
76  // get information about current chamber
77  Identifier mdt1_ID = (*mdt1)->identify();
78  bool mdt1_isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(mdt1_ID);
79  bool mdt1_isEndcap = m_idHelperSvc->mdtIdHelper().isEndcap(mdt1_ID);
80  int sector = m_idHelperSvc->sector(mdt1_ID);
81  int maxLayer = m_idHelperSvc->mdtIdHelper().tubeLayerMax(mdt1_ID);
82  int ML = m_idHelperSvc->mdtIdHelper().multilayer(mdt1_ID);
83 
84  // loop on hits inside the chamber
85  for (; mdt1 != mdtEnd; ++mdt1) {
86  if (Amg::error((*mdt1)->localCovariance(), Trk::locR) < m_errorCutOff) {
87  ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt1_ID) << " with too small error "
88  << Amg::error((*mdt1)->localCovariance(), Trk::locR));
89  continue;
90  }
91 
92  int tl1 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt1_ID);
93  if (tl1 == maxLayer) break; // require hits in at least 2 layers
94 
95  // loop on second hits
96  std::vector<const Muon::MdtPrepData*>::const_iterator mdt2 = (mdt1 + 1);
97  if (mdt2 == mdtEnd) continue;
98  Identifier mdt2_ID = (*mdt2)->identify();
99  for (; mdt2 != mdtEnd; ++mdt2) {
100  if (Amg::error((*mdt2)->localCovariance(), Trk::locR) < m_errorCutOff) {
101  ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt2_ID)
102  << " with too small error " << Amg::error((*mdt2)->localCovariance(), Trk::locR));
103  continue;
104  }
105 
106  // reject the bad tube combinations
107  int tl2 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt2_ID);
108  if (mdt1 == mdt2 || (tl2 - tl1) > 1 || (tl2 - tl1) < 0) continue;
109  if ((tl2 - tl1) == 0 && (m_idHelperSvc->mdtIdHelper().tube(mdt2_ID) -
110  m_idHelperSvc->mdtIdHelper().tube(mdt1_ID)) < 0) continue;
111  // reject bad hit separations
112  if (mdt1_isBarrel && std::abs((*mdt1)->globalPosition().z() - (*mdt2)->globalPosition().z()) > m_d12_max) continue;
113  if (mdt1_isEndcap && std::abs((*mdt1)->globalPosition().perp() - (*mdt2)->globalPosition().perp()) > m_d12_max) continue;
114 
115  // loop on third hits
116  std::vector<const Muon::MdtPrepData*>::const_iterator mdt3 = (mdt2 + 1);
117  if (mdt3 == mdtEnd) continue;
118  Identifier mdt3_ID = (*mdt3)->identify();
119  for (; mdt3 != mdtEnd; ++mdt3) {
120  if (Amg::error((*mdt3)->localCovariance(), Trk::locR) < m_errorCutOff) {
121  ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt3_ID)
122  << " with too small error " << Amg::error((*mdt3)->localCovariance(), Trk::locR));
123  continue;
124  }
125 
126  // reject the bad tube combinations
127  if (mdt1 == mdt3 || mdt2 == mdt3) continue;
128  int tl3 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt3_ID);
129  if ((tl3 - tl2) > 1 || (tl3 - tl2) < 0 || (tl3 - tl1) <= 0) continue;
130  if ((tl3 - tl2) == 0 && (m_idHelperSvc->mdtIdHelper().tube(mdt3_ID) -
131  m_idHelperSvc->mdtIdHelper().tube(mdt2_ID)) < 0) continue;
132  // reject bad hit separations
133  if (mdt1_isBarrel && std::abs((*mdt1)->globalPosition().z() - (*mdt3)->globalPosition().z()) > m_d13_max) continue;
134  if (mdt1_isEndcap && std::abs((*mdt1)->globalPosition().perp() - (*mdt3)->globalPosition().perp()) > m_d13_max) continue;
135 
136  // store and fit the good combinations
137  std::vector<const Muon::MdtPrepData*> mdts;
138  mdts.push_back((*mdt1));
139  mdts.push_back((*mdt2));
140  mdts.push_back((*mdt3));
141  std::vector<TrackletSegment> tmpSegs = TrackletSegmentFitter(mdts);
142  for (const TrackletSegment &tmpSeg : tmpSegs) mlsegments.push_back(tmpSeg);
143  } // end loop on mdt3
144  } // end loop on mdt2
145  } // end loop on mdt1
146 
147  // store the reconstructed segments according to station, ML and sector
148  // MS region decoded in MuonIdHelpers/MuonIdHelper.h
149  int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(mdt1_ID);
150  if (mdt1_isBarrel){
151  if (stationRegion == 0)
152  for (const TrackletSegment &mlsegment : mlsegments) segs[0][ML - 1][sector - 1].push_back(mlsegment); // barrel inner
153  else if (stationRegion == 2)
154  for (const TrackletSegment &mlsegment : mlsegments) segs[1][ML - 1][sector - 1].push_back(mlsegment); // barrel middle
155  else if (stationRegion == 3)
156  for (const TrackletSegment &mlsegment : mlsegments) segs[2][ML - 1][sector - 1].push_back(mlsegment); // barrel outer
157  }
158  else if (mdt1_isEndcap){
159  if (stationRegion == 0)
160  for (const TrackletSegment &mlsegment : mlsegments) segs[3][ML - 1][sector - 1].push_back(mlsegment); // endcap inner
161  else if (stationRegion == 2)
162  for (const TrackletSegment &mlsegment : mlsegments) segs[4][ML - 1][sector - 1].push_back(mlsegment); // endcap middle
163  else if (stationRegion == 3)
164  for (const TrackletSegment &mlsegment : mlsegments) segs[5][ML - 1][sector - 1].push_back(mlsegment); // endcap outer
165  }
166  else
167  ATH_MSG_WARNING("Found segments belonging to chamber " << m_idHelperSvc->mdtIdHelper().stationNameString(m_idHelperSvc->mdtIdHelper().stationName(mdt1_ID)) << " that have not been stored");
168  } // end loop on mdt chambers
169 
170  // Combine/remove duplicate segments
171  std::vector<TrackletSegment> CleanSegs[6][2][16];
172  for (int st = 0; st < 6; ++st) {
173  for (int ml = 0; ml < 2; ++ml) {
174  for (int sector = 0; sector < 16; ++sector) {
175  if (!segs[st][ml][sector].empty()) {
176  CleanSegs[st][ml][sector] = CleanSegments(segs[st][ml][sector]);
177  }
178  }
179  }
180  }
181 
182  // loop over TrackletSegments in barrel inner, middle, outer and endcap inner, middle, outer stations
183  for (int st = 0; st < 6; ++st) {
184  double DeltaAlphaCut = m_BarrelDeltaAlphaCut;
185  for (int sector = 0; sector < 16; ++sector) {
186  for (const TrackletSegment &ML1seg : CleanSegs[st][0][sector]) {
187  // Set the delta alpha cut depending on station type
188  const Identifier trkID = ML1seg.getIdentifier();
189  bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
190  bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
191  int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
192 
193  if (isBarrel){
194  if (stationRegion == 0){
195  if (isSmall) DeltaAlphaCut = m_BarrelDeltaAlphaCut; // default value for BIS
196  else DeltaAlphaCut = c_BIL / 750.0;
197  }
198  else if (stationRegion == 2){
199  if (isSmall) DeltaAlphaCut = c_BMS / 750.0;
200  else DeltaAlphaCut = c_BML / 750.0;
201  }
202  else if (stationRegion == 3){
203  if (isSmall) DeltaAlphaCut = m_BarrelDeltaAlphaCut; // default value for BOS
204  else DeltaAlphaCut = c_BOL / 750.0;
205  }
206  }
207  else{
208  DeltaAlphaCut = m_EndcapDeltaAlphaCut;
209  }
210 
211  // loop on ML2 segments from same sector
212  for (const TrackletSegment &ML2seg : CleanSegs[st][1][sector]) {
213  if (ML1seg.mdtChamber() != ML2seg.mdtChamber() || ML1seg.mdtChEta() != ML2seg.mdtChEta()) continue;
214 
215  double deltaAlpha = ML1seg.alpha() - ML2seg.alpha();
216  bool goodDeltab = DeltabCalc(ML1seg, ML2seg);
217  // select the good combinations
218  if (std::abs(deltaAlpha) < DeltaAlphaCut && goodDeltab) {
219  if (isBarrel) {
220  // barrel chambers
221  double charge_discriminant = deltaAlpha * ML1seg.globalPosition().z() * std::tan(ML1seg.alpha());
222  double charge = charge_discriminant < 0 ? -1 : 1;
223 
224  double pTot = TrackMomentum(ML1seg.getIdentifier(), deltaAlpha);
225  if (pTot < m_minpTot) continue;
226  if (pTot > m_maxpTot) {
227  // if we find a straight track, try to do a global refit to minimize the number of duplicates
228  charge = 0;
229  std::vector<const Muon::MdtPrepData*> mdts = ML1seg.mdtHitsOnTrack();
230  std::vector<const Muon::MdtPrepData*> mdts2 = ML2seg.mdtHitsOnTrack();
231  for (const Muon::MdtPrepData *mdt2 : mdts2) mdts.push_back(mdt2);
232  std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
233 
234  if (!CombinedSeg.empty()) {
235  // calculate momentum components & uncertainty
236  double Trk1overPErr = TrackMomentumError(CombinedSeg[0]);
237  double pT = pTot * std::sin(CombinedSeg[0].alpha());
238  double pz = pTot * std::cos(CombinedSeg[0].alpha());
239  Amg::Vector3D momentum(pT * std::cos(CombinedSeg[0].globalPosition().phi()),
240  pT * std::sin(CombinedSeg[0].globalPosition().phi()),
241  pz);
242  // create the error matrix
243  AmgSymMatrix(5) matrix;
244  matrix.setIdentity();
245  matrix(0, 0) = std::pow(CombinedSeg[0].rError(),2); // delta locR
246  matrix(1, 1) = std::pow(CombinedSeg[0].zError(),2); // delta locz
247  matrix(2, 2) = std::pow(0.00000000001,2); // delta phi (~0 because we explicitly rotate all tracklets into
248  // the middle of the chamber)
249  matrix(3, 3) = std::pow(CombinedSeg[0].alphaError(),2); // delta theta
250  matrix(4, 4) = std::pow(Trk1overPErr,2); // delta 1/p
251  Tracklet tmpTrk(CombinedSeg[0], momentum, matrix, charge);
252  ATH_MSG_DEBUG("Track " << tracklets.size() << " found with p = (" << momentum.x() << ", "
253  << momentum.y() << ", " << momentum.z()
254  << ") and |p| = " << tmpTrk.momentum().mag() << " MeV");
255  tracklets.push_back(tmpTrk);
256  }
257  } else {
258  // tracklet has a measurable momentum
259  double Trk1overPErr = TrackMomentumError(ML1seg, ML2seg);
260  double pT = pTot * std::sin(ML1seg.alpha());
261  double pz = pTot * std::cos(ML1seg.alpha());
262  Amg::Vector3D momentum(pT * std::cos(ML1seg.globalPosition().phi()),
263  pT * std::sin(ML1seg.globalPosition().phi()),
264  pz);
265  // create the error matrix
266  AmgSymMatrix(5) matrix;
267  matrix.setIdentity();
268  matrix(0, 0) = std::pow(ML1seg.rError(),2); // delta locR
269  matrix(1, 1) = std::pow(ML1seg.zError(),2); // delta locz
270  matrix(2, 2) = std::pow(0.00000000001,2); // delta phi (~0 because we explicitly rotate all tracks into the
271  // middle of the chamber)
272  matrix(3, 3) = std::pow(ML1seg.alphaError(),2); // delta theta
273  matrix(4, 4) = std::pow(Trk1overPErr,2); // delta 1/p
274  Tracklet tmpTrk(ML1seg, ML2seg, momentum, matrix, charge);
275  ATH_MSG_DEBUG("Track " << tracklets.size() << " found with p = (" << momentum.x() << ", "
276  << momentum.y() << ", " << momentum.z()
277  << ") and |p| = " << tmpTrk.momentum().mag() << " MeV");
278  tracklets.push_back(tmpTrk);
279  }
280  } // end barrel chamber selection
281  else if (!isBarrel) {
282  // endcap tracklets
283  // always straight tracklets (no momentum measurement possible)
284  std::vector<const Muon::MdtPrepData*> mdts = ML1seg.mdtHitsOnTrack();
285  std::vector<const Muon::MdtPrepData*> mdts2 = ML2seg.mdtHitsOnTrack();
286  for (const Muon::MdtPrepData *mdt2 : mdts2) mdts.push_back(mdt2);
287  std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
288 
289  if (!CombinedSeg.empty()) {
290  double charge = 0;
291  double pTot = m_straightTrackletpTot;
292  double pT = pTot * std::sin(CombinedSeg[0].alpha());
293  double pz = pTot * std::cos(CombinedSeg[0].alpha());
294  Amg::Vector3D momentum(pT * std::cos(CombinedSeg[0].globalPosition().phi()),
295  pT * std::sin(CombinedSeg[0].globalPosition().phi()),
296  pz);
297  // create the error matrix
298  AmgSymMatrix(5) matrix;
299  matrix.setIdentity();
300  matrix(0, 0) = std::pow(CombinedSeg[0].rError(),2); // delta locR
301  matrix(1, 1) = std::pow(CombinedSeg[0].zError(),2); // delta locz
302  matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle
303  // of the chamber)
304  matrix(3, 3) = std::pow(CombinedSeg[0].alphaError(),2); // delta theta
305  matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
306 
307  Tracklet tmpTrk(CombinedSeg[0], momentum, matrix, charge);
308  tracklets.push_back(tmpTrk);
309  }
310  } // end endcap tracklet selection
311 
312  } // end tracklet selection (delta alpha & delta b)
313 
314  } // end loop on ML2 segments
315  } // end loop on ML1 segments
316  } // end loop on sectors
317  } // end loop on stations
318 
319  // Resolve any ambiguous tracklets
320  tracklets = ResolveAmbiguousTracklets(tracklets);
321 
322  // convert from tracklets to Trk::Tracks
323  convertToTrackParticles(tracklets, container);
324 
325  return StatusCode::SUCCESS;
326  }
327 
328  //** ----------------------------------------------------------------------------------------------------------------- **//
329 
330  void MSVertexTrackletTool::convertToTrackParticles(std::vector<Tracklet>& tracklets,
331  SG::WriteHandle<xAOD::TrackParticleContainer>& container) {
332  // convert tracklets to xAOD::TrackParticle and store in a TrackCollection
333  for (Tracklet &tracklet : tracklets) {
334  xAOD::TrackParticle* trackparticle = new xAOD::TrackParticle();
335  tracklet.setTrackParticle(trackparticle);
336  container->push_back(trackparticle);
337 
338  AmgSymMatrix(5) covariance{tracklet.errorMatrix()};
339  auto MyPerigee(std::make_unique<Trk::Perigee>(tracklet.globalPosition(), tracklet.momentum(), tracklet.charge(), Trk::PerigeeSurface(Amg::Vector3D::Zero()), covariance));
340 
341  // fill the xAOD::TrackParticle with the tracklet content
342  trackparticle->setDefiningParameters(MyPerigee->parameters()[Trk::d0], MyPerigee->parameters()[Trk::z0],
343  MyPerigee->parameters()[Trk::phi0], MyPerigee->parameters()[Trk::theta],
344  MyPerigee->parameters()[Trk::qOverP]);
345  trackparticle->setFitQuality(1., (float)tracklet.mdtHitsOnTrack().size());
346  trackparticle->setTrackProperties(xAOD::TrackProperties::LowPtTrack);
347  std::vector<float> covMatrixVec;
348  Amg::compress(covariance, covMatrixVec);
349  trackparticle->setDefiningParametersCovMatrixVec(covMatrixVec);
350  }
351  return;
352  }
353 
354  //** ----------------------------------------------------------------------------------------------------------------- **//
355 
356  bool MSVertexTrackletTool::IgnoreMDTChamber(const Muon::MdtPrepData* mdtHit) const {
357  // return true if the MDT hit is in a chamber to be ignored. These hits are then not used to reconstruct tracklets.
358 
359  bool ignore = false;
360  int stName = m_idHelperSvc->mdtIdHelper().stationName(mdtHit->identify());
361  int stEta = m_idHelperSvc->mdtIdHelper().stationEta(mdtHit->identify());
362 
363  // Doesn't consider hits belonging to chambers BEE, EEL and EES
364  if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BEE") ||
365  stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("EEL") ||
366  stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("EES")) ignore = true;
367 
368  // Doesn't consider hits belonging to chambers BIS7/8
369  if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BIS") && std::abs(stEta) >= 7) ignore = true;
370 
371  // Doesn't consider hits belonging to BME or BMG chambers
372  if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BME") ||
373  stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BMG")) ignore = true;
374 
375  return ignore;
376  }
377 
378  //** ----------------------------------------------------------------------------------------------------------------- **//
379 
380  int MSVertexTrackletTool::SortMDThits(std::vector<std::vector<const Muon::MdtPrepData*> >& SortedMdt, const EventContext& ctx) const {
381  SortedMdt.clear();
382  int nMDT(0);
383 
384  SG::ReadHandle<Muon::MdtPrepDataContainer> mdtTES(m_mdtTESKey, ctx);
385  if (!mdtTES.isValid()) {
386  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Muon::MdtPrepDataContainer with key MDT_DriftCircles was not retrieved" << endmsg;
387  return 0;
388  } else {
389  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Muon::MdtPrepDataContainer with key MDT_DriftCircles retrieved" << endmsg;
390  }
391 
392  // iterators over collections, a collection corresponds to a chamber
393  for (const Muon::MdtPrepDataCollection* MDTch : *mdtTES){
394  if (MDTch->empty()) continue;
395  if (IgnoreMDTChamber(*(MDTch->begin()))) continue;
396 
397  // sort per multi layer
398  std::vector<const Muon::MdtPrepData*> hitsML1;
399  std::vector<const Muon::MdtPrepData*> hitsML2;
400 
401  // loop on mdt hits in the current chamber
402  for (const Muon::MdtPrepData* mdt : *MDTch) {
403  // Removes noisy hits
404  if (mdt->adc() < 50) continue;
405  // Removes dead modules or out of time hits
406  if (mdt->status() != Muon::MdtStatusDriftTime) continue;
407  // Removes tubes out of readout during drift time or with unphysical errors
408  if (mdt->localPosition()[Trk::locR] == 0.) continue;
409  if (mdt->localCovariance()(Trk::locR, Trk::locR) < 1e-6) {
410  ATH_MSG_WARNING("Found MDT with unphysical error " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt->identify())
411  << " cov " << mdt->localCovariance()(Trk::locR, Trk::locR));
412  continue;
413  }
414  ++nMDT;
415 
416  // sort per multi layer
417  if (m_idHelperSvc->mdtIdHelper().multilayer(mdt->identify()) == 1)
418  hitsML1.push_back(mdt);
419  else
420  hitsML2.push_back(mdt);
421 
422  } // end MdtPrepDataCollection
423 
424  // add
425  addMDTHits(hitsML1, SortedMdt);
426  addMDTHits(hitsML2, SortedMdt);
427  } // end MdtPrepDataContainer
428 
429  return nMDT;
430  }
431 
432  void MSVertexTrackletTool::addMDTHits(std::vector<const Muon::MdtPrepData*>& hits,
433  std::vector<std::vector<const Muon::MdtPrepData*> >& SortedMdt) const {
434  if (hits.empty()) return;
435 
436  // calculate number of hits in ML
437  int ntubes = hits.front()->detectorElement()->getNLayers() * hits.front()->detectorElement()->getNtubesperlayer();
438  if (hits.size() > 0.75 * ntubes) return;
439  std::sort(hits.begin(), hits.end(), [this](const Muon::MdtPrepData* mprd1, const Muon::MdtPrepData* mprd2) -> bool {
440  if (m_idHelperSvc->mdtIdHelper().tubeLayer(mprd1->identify()) > m_idHelperSvc->mdtIdHelper().tubeLayer(mprd2->identify()))
441  return false;
442  if (m_idHelperSvc->mdtIdHelper().tubeLayer(mprd1->identify()) < m_idHelperSvc->mdtIdHelper().tubeLayer(mprd2->identify()))
443  return true;
444  if (m_idHelperSvc->mdtIdHelper().tube(mprd1->identify()) < m_idHelperSvc->mdtIdHelper().tube(mprd2->identify())) return true;
445  return false;
446  }); // sort the MDTs by layer and tube number
447 
448  SortedMdt.push_back(hits);
449  }
450 
451  //** ----------------------------------------------------------------------------------------------------------------- **//
452 
453  std::vector<TrackletSegment> MSVertexTrackletTool::TrackletSegmentFitter(const std::vector<const Muon::MdtPrepData*>& mdts) const {
454  // fits TrackletSegments from an as compatible identified set of MDT hits
455  // create the segment seeds
456  std::vector<std::pair<double, double> > SeedParams = SegSeeds(mdts);
457  // fit the segments
458  std::vector<TrackletSegment> segs = TrackletSegmentFitterCore(mdts, SeedParams);
459 
460  return segs;
461  }
462 
463  //** ----------------------------------------------------------------------------------------------------------------- **//
464 
465  std::vector<std::pair<double, double> > MSVertexTrackletTool::SegSeeds(const std::vector<const Muon::MdtPrepData*>& mdts) const {
466  std::vector<std::pair<double, double> > SeedParams;
467  // create seeds by drawing the 4 possible lines tangent to the two outermost drift circles
468  // see http://cds.cern.ch/record/620198 (section 4.3) for description of the algorithm
469  // keep all seeds which satisfy the criterion: residual(mdt 2) < m_SeedResidual
470  // NOTE: here there is an assumption that each MDT has a radius of 30mm
471  // -- needs to be revisited when the small tubes in sectors 12 & 14 are installed
472  double x1 = mdts.front()->globalPosition().z();
473  double y1 = mdts.front()->globalPosition().perp();
474  double r1 = std::abs(mdts.front()->localPosition()[Trk::locR]);
475 
476  double x2 = mdts.back()->globalPosition().z();
477  double y2 = mdts.back()->globalPosition().perp();
478  double r2 = std::abs(mdts.back()->localPosition()[Trk::locR]);
479 
480  double DeltaX = x2 - x1;
481  double DeltaY = y2 - y1;
482  double DistanceOfCenters = std::hypot(DeltaX, DeltaY);
483  if (DistanceOfCenters < 30) return SeedParams;
484  double Alpha0 = std::acos(DeltaX / DistanceOfCenters);
485 
486  // First seed
487  double phi = mdts.front()->globalPosition().phi();
488  double RSum = r1 + r2;
489  if (RSum > DistanceOfCenters) return SeedParams;
490  double Alpha1 = std::asin(RSum / DistanceOfCenters);
491  double line_theta = Alpha0 + Alpha1;
492  double z_line = x1 + r1 * std::sin(line_theta);
493  double rho_line = y1 - r1 * std::cos(line_theta);
494 
495  Amg::Vector3D gPos1(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
496  Amg::Vector3D gDir(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
497  Amg::Vector3D globalDir1(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
498  double gSlope1 = (globalDir1.perp() / globalDir1.z());
499  double gInter1 = gPos1.perp() - gSlope1 * gPos1.z();
500  double resid = SeedResiduals(mdts, gSlope1, gInter1);
501  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope1, gInter1);
502  // Second seed
503  line_theta = Alpha0 - Alpha1;
504  z_line = x1 - r1 * std::sin(line_theta);
505  rho_line = y1 + r1 * std::cos(line_theta);
506  Amg::Vector3D gPos2(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
507  Amg::Vector3D globalDir2(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
508  double gSlope2 = (globalDir2.perp() / globalDir2.z());
509  double gInter2 = gPos2.perp() - gSlope2 * gPos2.z();
510  resid = SeedResiduals(mdts, gSlope2, gInter2);
511  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope2, gInter2);
512 
513  double Alpha2 = std::asin(std::abs(r2 - r1) / DistanceOfCenters);
514  if (r1 < r2) {
515  // Third seed
516  line_theta = Alpha0 + Alpha2;
517  z_line = x1 - r1 * std::sin(line_theta);
518  rho_line = y1 + r1 * std::cos(line_theta);
519 
520  Amg::Vector3D gPos3(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
521  Amg::Vector3D globalDir3(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
522  double gSlope3 = (globalDir3.perp() / globalDir3.z());
523  double gInter3 = gPos3.perp() - gSlope3 * gPos3.z();
524  resid = SeedResiduals(mdts, gSlope3, gInter3);
525  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope3, gInter3);
526 
527  // Fourth seed
528  line_theta = Alpha0 - Alpha2;
529  z_line = x1 + r1 * std::sin(line_theta);
530  rho_line = y1 - r1 * std::cos(line_theta);
531 
532  Amg::Vector3D gPos4(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
533  Amg::Vector3D globalDir4(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
534  double gSlope4 = (globalDir4.perp() / globalDir4.z());
535  double gInter4 = gPos4.perp() - gSlope4 * gPos4.z();
536  resid = SeedResiduals(mdts, gSlope4, gInter4);
537  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope4, gInter4);
538  } else {
539  // Third seed
540  line_theta = Alpha0 + Alpha2;
541  z_line = x1 + r1 * std::sin(line_theta);
542  rho_line = y1 - r1 * std::cos(line_theta);
543 
544  Amg::Vector3D gPos3(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
545  Amg::Vector3D globalDir3(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
546  double gSlope3 = (globalDir3.perp() / globalDir3.z());
547  double gInter3 = gPos3.perp() - gSlope3 * gPos3.z();
548  resid = SeedResiduals(mdts, gSlope3, gInter3);
549  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope3, gInter3);
550 
551  // Fourth seed
552  line_theta = Alpha0 - Alpha2;
553  z_line = x1 - r1 * std::sin(line_theta);
554  rho_line = y1 + r1 * std::cos(line_theta);
555 
556  Amg::Vector3D gPos4(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
557  Amg::Vector3D globalDir4(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
558  double gSlope4 = (globalDir4.perp() / globalDir4.z());
559  double gInter4 = gPos4.perp() - gSlope4 * gPos4.z();
560  resid = SeedResiduals(mdts, gSlope4, gInter4);
561  if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope4, gInter4);
562  }
563 
564  return SeedParams;
565  }
566 
567  //** ----------------------------------------------------------------------------------------------------------------- **//
568 
569  double MSVertexTrackletTool::SeedResiduals(const std::vector<const Muon::MdtPrepData*>& mdts, double slope, double inter) {
570  // calculate the residual of the MDTs not used to create the seed
571  double resid = 0;
572  for (const Muon::MdtPrepData *mdt : mdts) {
573  double mdtR = mdt->globalPosition().perp();
574  double mdtZ = mdt->globalPosition().z();
575  double res = std::abs((mdt->localPosition()[Trk::locR] - std::abs((mdtR - inter - slope * mdtZ) / std::hypot(slope, 1))) /
576  (Amg::error(mdt->localCovariance(), Trk::locR)));
577  if (res > resid) resid = res;
578  }
579  return resid;
580  }
581 
582  //** ----------------------------------------------------------------------------------------------------------------- **//
583 
584  std::vector<TrackletSegment> MSVertexTrackletTool::TrackletSegmentFitterCore(const std::vector<const Muon::MdtPrepData*>& mdts,
585  const std::vector<std::pair<double, double> >& SeedParams) const {
586  std::vector<TrackletSegment> segs;
587 
588  Identifier mdtID = mdts.at(0)->identify();
589  for (const std::pair<double,double> &SeedParam : SeedParams) {
590  // Min chi^2 fit from "Precision of the ATLAS Muon Spectrometer" -- M. Woudstra
591  // http://cds.cern.ch/record/620198?ln=en (section 4.3)
592  double chi2(0);
593  double s(0), sz(0), sy(0);
594  // loop on the mdt hits, find the weighted center
595  for (const Muon::MdtPrepData *prd : mdts) {
596  // Tell clang to optimize assuming that FP exceptions can trap.
597  // Otherwise, it can vectorize the division, which can lead to
598  // spurious division-by-zero traps from unused vector lanes.
599  CXXUTILS_TRAPPING_FP;
600  const double mdt_y = std::hypot(prd->globalPosition().x(), prd->globalPosition().y());
601  const double mdt_z = prd->globalPosition().z();
602  const double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
603  s += 1 / sigma2;
604  sz += mdt_z / sigma2;
605  sy += mdt_y / sigma2;
606  }
607  const double yc = sy / s;
608  const double zc = sz / s;
609 
610  // Find the initial parameters of the fit
611  double alpha = std::atan2(SeedParam.first, 1.0);
612  if (alpha < 0) alpha += M_PI;
613  double dalpha = 0;
614  double d = (SeedParam.second - yc + zc * SeedParam.first) * std::cos(alpha);
615  double dd = 0;
616 
617  // require segments to point to the second ML
618  if (std::abs(std::cos(alpha)) > 0.97 && (m_idHelperSvc->mdtIdHelper().isBarrel(mdtID))) continue;
619  if (std::abs(std::cos(alpha)) < 0.03 && (m_idHelperSvc->mdtIdHelper().isEndcap(mdtID))) continue;
620 
621  // calculate constants used in the fit
622  double sPyy(0), sPyz(0), sPyyzz(0);
623  for (const Muon::MdtPrepData *prd : mdts) {
624  double mdt_y = std::hypot(prd->globalPosition().x(), prd->globalPosition().y());
625  double mdt_z = prd->globalPosition().z();
626  double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
627  sPyy += std::pow(mdt_y-yc,2) / sigma2;
628  sPyz += (mdt_y - yc) * (mdt_z - zc) / sigma2;
629  sPyyzz += ((mdt_y - yc) - (mdt_z - zc)) * ((mdt_y - yc) + (mdt_z - zc)) / sigma2;
630  }
631 
632  // iterative fit
633  int Nitr = 0;
634  double deltaAlpha = 0;
635  double deltad = 0;
636  while (true) {
637  double sumRyi(0), sumRzi(0), sumRi(0);
638  chi2 = 0;
639  ++Nitr;
640  const double cos_a = std::cos(alpha);
641  const double sin_a = std::sin(alpha);
642  for (const Muon::MdtPrepData *prd : mdts) {
643  double mdt_y = prd->globalPosition().perp();
644  double mdt_z = prd->globalPosition().z();
645  double yPi = -(mdt_z - zc) * sin_a + (mdt_y - yc) * cos_a - d;
646  double signR = yPi >= 0 ? -1. : 1;
647  double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
648  double ri = signR * prd->localPosition()[Trk::locR];
649  sumRyi += ri * (mdt_y - yc) / sigma2;
650  sumRzi += ri * (mdt_z - zc) / sigma2;
651  sumRi += ri / sigma2;
652  chi2 += std::pow(yPi+ri,2) / sigma2;
653  }
654  double bAlpha = -1 * sPyz + cos_a * (sin_a * sPyyzz + 2 * cos_a * sPyz + sumRzi) + sin_a * sumRyi;
655  double AThTh = sPyy + cos_a * (2 * sin_a * sPyz - cos_a * sPyyzz);
659  if (std::abs(AThTh) < 1.e-7) break;
660  // the new alpha & d parameters
661  double alphaNew = alpha + bAlpha / AThTh;
662  double dNew = sumRi / s;
663  // the errors
664  dalpha = std::sqrt(1 / std::abs(AThTh));
665  dd = std::sqrt(1 / s);
666  deltaAlpha = std::abs(alphaNew - alpha);
667  deltad = std::abs(d - dNew);
668  // test if the new segment is different than the previous
669  if (deltaAlpha < 5.e-7 && deltad < 5.e-6) break;
670  alpha = alphaNew;
671  d = dNew;
672  // Guard against infinite loops
673  if (Nitr > 10) break;
674  } // end while loop
675 
676  // find the chi^2 probability of the segment
677  double chi2Prob = TMath::Prob(chi2, mdts.size() - 2);
678  // keep only "good" segments
679  if (chi2Prob > m_minSegFinderChi2) {
680  double z0 = zc - d * std::sin(alpha);
681  double dz0 = std::hypot(dd*std::sin(alpha), d*dalpha*std::cos(alpha));
682  double y0 = yc + d * std::cos(alpha);
683  double dy0 = std::hypot(dd*std::cos(alpha), d*dalpha*std::sin(alpha));
684  // find the hit pattern, which side of the wire did the particle pass? (1==Left, 2==Right)
685  /*
686  ( )(/O)( )
687  (./)( )( ) == RRL == 221
688  (O/)( )( )
689  */
690  int pattern(0);
691  if (mdts.size() > 8)
692  pattern = -1; // with more then 8 MDTs the pattern is unique
693  else {
694  for (unsigned int k = 0; k < mdts.size(); ++k) {
695  int base = std::pow(10, k);
696  double mdtR = std::hypot(mdts.at(k)->globalPosition().x(), mdts.at(k)->globalPosition().y());
697  double mdtZ = mdts.at(k)->globalPosition().z();
698  double zTest = (mdtR - y0) / std::tan(alpha) + z0 - mdtZ;
699  if (zTest > 0)
700  pattern += 2 * base;
701  else
702  pattern += base;
703  }
704  }
705 
706  // find the position of the tracklet in the global frame
707  double mdtPhi = mdts.at(0)->globalPosition().phi();
708  Amg::Vector3D segpos(y0 * std::cos(mdtPhi), y0 * std::sin(mdtPhi), z0);
709  // create the tracklet
710  TrackletSegment MyTrackletSegment{m_idHelperSvc.get(), mdts, segpos, alpha, dalpha, dy0, dz0, pattern};
711  segs.push_back(MyTrackletSegment);
712  if (pattern == -1) break; // stop if we find a segment with more than 8 hits (guaranteed to be unique!)
713  }
714  } // end loop on segment seeds
715 
716  // in case more than 1 segment is reconstructed, check if there are duplicates using the hit patterns
717  if (segs.size() > 1) {
718  std::vector<TrackletSegment> tmpSegs;
719  for (unsigned int i1 = 0; i1 < segs.size(); ++i1) {
720  bool isUnique = true;
721  int pattern1 = segs.at(i1).getHitPattern();
722  for (unsigned int i2 = (i1 + 1); i2 < segs.size(); ++i2) {
723  if (pattern1 == -1) break;
724  int pattern2 = segs.at(i2).getHitPattern();
725  if (pattern1 == pattern2) isUnique = false;
726  }
727  if (isUnique) tmpSegs.push_back(segs.at(i1));
728  }
729  segs = tmpSegs;
730  }
731 
732  // return the unique segments
733  return segs;
734  }
735 
736  //** ----------------------------------------------------------------------------------------------------------------- **//
737 
738  std::vector<TrackletSegment> MSVertexTrackletTool::CleanSegments(const std::vector<TrackletSegment>& Segs) const {
739  std::vector<TrackletSegment> CleanSegs;
740  std::vector<TrackletSegment> segs = Segs; // set of segments to perform cleaning on
741  bool keepCleaning(true);
742  int nItr(0);
743 
744  while (keepCleaning) {
745  ++nItr;
746  keepCleaning = false;
747 
748  for (std::vector<TrackletSegment>::iterator it = segs.begin(); it != segs.end(); ++it) {
749  if (it->isCombined()) continue;
750  std::vector<TrackletSegment> segsToCombine;
751  double tanTh1 = std::tan(it->alpha());
752  double r1 = it->globalPosition().perp();
753  double zi1 = it->globalPosition().z() - r1 / tanTh1;
754  // find all segments with similar parameters & attempt to combine
755  for (std::vector<TrackletSegment>::iterator sit = (it + 1); sit != segs.end(); ++sit) {
756  if (sit->isCombined()) continue;
757  if (it->mdtChamber() != sit->mdtChamber()) continue; // require the segments are in the same chamber
758  if ((it->mdtChEta()) * (sit->mdtChEta()) < 0) continue; // check both segments are on the same side of the detector
759  if (it->mdtChPhi() != sit->mdtChPhi()) continue; // in the same sector
760  if (std::abs(it->alpha() - sit->alpha()) > 0.005) continue; // same trajectory
761  double tanTh2 = std::tan(sit->alpha());
762  double r2 = sit->globalPosition().perp();
763  double zi2 = sit->globalPosition().z() - r2 / tanTh2;
764  // find the distance at the midpoint between the two segments
765  double rmid = (r1 + r2) / 2.;
766  double z1 = rmid / tanTh1 + zi1;
767  double z2 = rmid / tanTh2 + zi2;
768  double zdist = std::abs(z1 - z2);
769  if (zdist < 0.5) {
770  segsToCombine.push_back(*sit);
771  sit->isCombined(true);
772  }
773  } // end sit loop
774 
775  // if the segment is unique, keep it
776  if (segsToCombine.empty()) {
777  CleanSegs.push_back(*it);
778  }
779  // else, combine all like segments & refit
780  else if (!segsToCombine.empty()) {
781  // create a vector of all unique MDT hits in the segments
782  std::vector<const Muon::MdtPrepData*> mdts = it->mdtHitsOnTrack();
783  for (const TrackletSegment &seg : segsToCombine) {
784  std::vector<const Muon::MdtPrepData*> tmpmdts = seg.mdtHitsOnTrack();
785  for (const Muon::MdtPrepData *tmpprd : tmpmdts){
786  bool isNewHit(true);
787  for (const Muon::MdtPrepData *tmpprd2 : mdts){
788  if (tmpprd->identify() == tmpprd2->identify()) {
789  isNewHit = false;
790  break;
791  }
792  }
793  if (isNewHit && Amg::error(tmpprd->localCovariance(), Trk::locR) > m_errorCutOff) mdts.push_back(tmpprd);
794  }
795  } // end segsToCombine loop
796 
797  // only need to combine if there are extra hits added to the first segment
798  if (mdts.size() > it->mdtHitsOnTrack().size()) {
799  std::vector<TrackletSegment> refitsegs = TrackletSegmentFitter(mdts);
800  // if the refit fails, what to do?
801  if (refitsegs.empty()) {
802  if (segsToCombine.size() == 1) {
803  segsToCombine[0].isCombined(false);
804  CleanSegs.push_back(*it);
805  CleanSegs.push_back(segsToCombine[0]);
806  } else {
807  // loop on the mdts and count the number of segments that share that hit
808  std::vector<int> nSeg;
809  for (unsigned int i = 0; i < mdts.size(); ++i) {
810  nSeg.push_back(0);
811  // hit belongs to the first segment
812  for (unsigned int k = 0; k < it->mdtHitsOnTrack().size(); ++k) {
813  if (it->mdtHitsOnTrack()[k]->identify() == mdts[i]->identify()) {
814  ++nSeg[i];
815  break;
816  }
817  }
818  // hit belongs to one of the duplicate segments
819  for (unsigned int k = 0; k < segsToCombine.size(); ++k) {
820  for (unsigned int m = 0; m < segsToCombine[k].mdtHitsOnTrack().size(); ++m) {
821  if (segsToCombine[k].mdtHitsOnTrack()[m]->identify() == mdts[i]->identify()) {
822  ++nSeg[i];
823  break;
824  }
825  } // end loop on mdtHitsOnTrack
826  } // end loop on segsToCombine
827  } // end loop on mdts
828 
829  // loop over the duplicates and remove the MDT used by the fewest segments until the fit converges
830  bool keeprefitting(true);
831  int nItr2(0);
832  while (keeprefitting) {
833  ++nItr2;
834  int nMinSeg(nSeg[0]);
835  const Muon::MdtPrepData* minmdt = mdts[0];
836  std::vector<int> nrfsegs;
837  std::vector<const Muon::MdtPrepData*> refitmdts;
838  // loop on MDTs, identify the overlapping set of hits
839  for (unsigned int i = 1; i < mdts.size(); ++i) {
840  if (nSeg[i] < nMinSeg) {
841  refitmdts.push_back(minmdt);
842  nrfsegs.push_back(nMinSeg);
843  minmdt = mdts[i];
844  nMinSeg = nSeg[i];
845  } else {
846  refitmdts.push_back(mdts[i]);
847  nrfsegs.push_back(nSeg[i]);
848  }
849  }
850  // reset the list of MDTs & the minimum number of segments an MDT must belong to
851  mdts = refitmdts;
852  nSeg = nrfsegs;
853  // try to fit the new set of MDTs
854  refitsegs = TrackletSegmentFitter(mdts);
855  if (!refitsegs.empty()) {
856  for (const TrackletSegment &refitseg : refitsegs) CleanSegs.push_back(refitseg);
857  keeprefitting = false; // stop refitting if segments are found
858  } else if (mdts.size() <= 3) {
859  CleanSegs.push_back(*it);
860  keeprefitting = false;
861  }
862  if (nItr2 > 10) break;
863  } // end while
864  }
865  } else {
866  keepCleaning = true;
867  for (const TrackletSegment &refitseg : refitsegs) CleanSegs.push_back(refitseg);
868  }
869  }
870  // if there are no extra MDT hits, keep only the first segment as unique
871  else
872  CleanSegs.push_back(*it);
873  }
874  } // end it loop
875  if (keepCleaning) {
876  segs = CleanSegs;
877  CleanSegs.clear();
878  }
879  if (nItr > 10) break;
880  } // end while
881 
882  return CleanSegs;
883  }
884 
885  //** ----------------------------------------------------------------------------------------------------------------- **//
886 
887  bool MSVertexTrackletTool::DeltabCalc(const TrackletSegment& ML1seg, const TrackletSegment& ML2seg) const {
888  double ChMid = (ML1seg.getChMidPoint() + ML2seg.getChMidPoint()) / 2.0;
889  // Calculate the Delta b (see http://inspirehep.net/record/1266438)
890  double mid1(100), mid2(1000);
891  double deltab(100);
892  if (m_idHelperSvc->mdtIdHelper().isBarrel(ML1seg.getIdentifier())) {
893  // delta b in the barrel
894  mid1 = (ChMid - ML1seg.globalPosition().perp()) / std::tan(ML1seg.alpha()) + ML1seg.globalPosition().z();
895  mid2 = (ChMid - ML2seg.globalPosition().perp()) / std::tan(ML2seg.alpha()) + ML2seg.globalPosition().z();
896  double r01 = ML1seg.globalPosition().perp() - ML1seg.globalPosition().z() * std::tan(ML1seg.alpha());
897  double r02 = ML2seg.globalPosition().perp() - ML2seg.globalPosition().z() * std::tan(ML2seg.alpha());
898  deltab = (mid2 * std::tan(ML1seg.alpha()) - ChMid + r01) / (std::hypot(1,std::tan(ML1seg.alpha())));
899  double deltab2 = (mid1 * std::tan(ML2seg.alpha()) - ChMid + r02) / (std::hypot(1,std::tan(ML2seg.alpha())));
900  if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
901  } else {
902  // delta b in the endcap
903  mid1 = ML1seg.globalPosition().perp() + std::tan(ML1seg.alpha()) * (ChMid - ML1seg.globalPosition().z());
904  mid2 = ML2seg.globalPosition().perp() + std::tan(ML2seg.alpha()) * (ChMid - ML2seg.globalPosition().z());
905  double z01 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
906  double z02 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
907  deltab = (mid2 / std::tan(ML1seg.alpha()) - ChMid + z01) / (std::hypot(1,1/std::tan(ML1seg.alpha())));
908  double deltab2 = (mid1 / std::tan(ML2seg.alpha()) - ChMid + z02) / (std::hypot(1,1/std::tan(ML2seg.alpha())));
909  if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
910  }
911 
912  // calculate the maximum allowed Delta b based on delta alpha uncertainties and ML spacing
913  double dbmax = 5 * std::abs(ChMid - ML1seg.getChMidPoint()) * std::hypot(ML1seg.alphaError(), ML2seg.alphaError());
914  if (dbmax > m_maxDeltabCut) dbmax = m_maxDeltabCut;
915  return std::abs(deltab) < dbmax;
916  }
917 
918  //** ----------------------------------------------------------------------------------------------------------------- **//
919 
920  double MSVertexTrackletTool::TrackMomentum(const Identifier trkID, const double deltaAlpha) const {
921  // p = k/delta_alpha
922  bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
923  bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
924  int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
925 
926  double dalpha = std::abs(deltaAlpha);
927  double pTot = m_straightTrackletpTot;
928  if (isBarrel){
929  if (stationRegion == 0){
930  if (isSmall) pTot = m_straightTrackletpTot; // default value for BIS
931  else pTot = c_BIL / dalpha;
932  }
933  else if (stationRegion == 2){
934  if (isSmall) pTot = c_BMS / dalpha;
935  else pTot = c_BML / dalpha;
936  }
937  else if (stationRegion == 3){
938  if (isSmall) pTot = m_straightTrackletpTot; // default value for BOS
939  else pTot = c_BOL / dalpha;
940  }
941  }
942 
943  if (pTot > m_maxpTot) pTot = m_straightTrackletpTot;
944 
945  return pTot;
946  }
947 
948  //** ----------------------------------------------------------------------------------------------------------------- **//
949 
950  double MSVertexTrackletTool::TrackMomentumError(const TrackletSegment& ml1, const TrackletSegment& ml2) const{
951  // uncertainty on 1/p
952  const Identifier trkID = ml1.getIdentifier();
953  bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
954  bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
955  int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
956 
957  double dalpha = std::hypot(ml1.alphaError(), ml2.alphaError());
958  double pErr = dalpha / c_BML;
959  if (isBarrel){
960  if (stationRegion == 0){
961  if (isSmall) pErr = dalpha / c_BML; // default value for BIS
962  else pErr = dalpha / c_BIL;
963  }
964  else if (stationRegion == 2){
965  if (isSmall) pErr = dalpha / c_BMS;
966  else pErr = dalpha / c_BML;
967  }
968  else if (stationRegion == 3){
969  if (isSmall) pErr = dalpha / c_BML; // default value for BOS
970  else pErr = dalpha / c_BOL;
971  }
972  }
973 
974  return pErr;
975  }
976 
977  //** ----------------------------------------------------------------------------------------------------------------- **//
978 
979  double MSVertexTrackletTool::TrackMomentumError(const TrackletSegment& ml1) const {
980  // uncertainty in 1/p
981  const Identifier trkID = ml1.getIdentifier();
982  bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
983  bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
984  int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
985 
986  double dalpha = std::abs(ml1.alphaError());
987  double pErr = dalpha / c_BML;
988  if (isBarrel){
989  if (stationRegion == 0){
990  if (isSmall) pErr = dalpha / c_BML; // default value for BIS
991  else pErr = dalpha / c_BIL;
992  }
993  else if (stationRegion == 2){
994  if (isSmall) pErr = dalpha / c_BMS;
995  else pErr = dalpha / c_BML;
996  }
997  else if (stationRegion == 3){
998  if (isSmall) pErr = dalpha / c_BML; // default value for BOS
999  else pErr = dalpha / c_BOL;
1000  }
1001  }
1002 
1003  return pErr;
1004  }
1005 
1006  //** ----------------------------------------------------------------------------------------------------------------- **//
1007 
1008  std::vector<Tracklet> MSVertexTrackletTool::ResolveAmbiguousTracklets(std::vector<Tracklet>& tracks) const {
1009  ATH_MSG_DEBUG("In ResolveAmbiguousTracks");
1010  // considering only tracklets with the number of associated hits
1011  // being more than 3/4 the number of layers in the MS chamber
1012 
1013  if (m_tightTrackletRequirement) {
1014  std::vector<Tracklet> myTracks = tracks;
1015  tracks.clear();
1016  for (const Tracklet &track : myTracks) {
1017  Identifier id1 = track.getML1seg().mdtHitsOnTrack().at(0)->identify();
1018  Identifier id2 = track.getML2seg().mdtHitsOnTrack().at(0)->identify();
1019  int nLayerML1 = m_idHelperSvc->mdtIdHelper().tubeLayerMax(id1);
1020  int nLayerML2 = m_idHelperSvc->mdtIdHelper().tubeLayerMax(id2);
1021  double ratio = (double)(track.mdtHitsOnTrack().size()) / (nLayerML1 + nLayerML2);
1022  if (ratio > 0.75) tracks.push_back(track);
1023  }
1024  }
1025 
1026  std::vector<Tracklet> UniqueTracks;
1027  std::vector<unsigned int> AmbigTrks; // indices of ambigious tracklets
1028  for (unsigned int tk1 = 0; tk1 < tracks.size(); ++tk1) {
1029  int nShared = 0;
1030  // check if any Ambiguity has been broken
1031  bool isResolved = false;
1032  for (unsigned int AmbigTrksIdx : AmbigTrks) {
1033  if (tk1 == AmbigTrksIdx) {
1034  isResolved = true;
1035  break;
1036  }
1037  }
1038  if (isResolved) continue;
1039  std::vector<Tracklet> AmbigTracks;
1040  AmbigTracks.push_back(tracks.at(tk1));
1041  // get a point on the track
1042  double Trk1ML1R = tracks.at(tk1).getML1seg().globalPosition().perp();
1043  double Trk1ML1Z = tracks.at(tk1).getML1seg().globalPosition().z();
1044  double Trk1ML2R = tracks.at(tk1).getML2seg().globalPosition().perp();
1045  double Trk1ML2Z = tracks.at(tk1).getML2seg().globalPosition().z();
1046 
1047  Identifier tk1ID = tracks.at(tk1).muonIdentifier();
1048  bool tk1_isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(tk1ID);
1049  bool tk1_isEndcap = m_idHelperSvc->mdtIdHelper().isEndcap(tk1ID);
1050 
1051  // loop over the rest of the tracks and find any amibuities
1052  for (unsigned int tk2 = (tk1 + 1); tk2 < tracks.size(); ++tk2) {
1053  if (tracks.at(tk1).mdtChamber() == tracks.at(tk2).mdtChamber() && tracks.at(tk1).mdtChPhi() == tracks.at(tk2).mdtChPhi() &&
1054  (tracks.at(tk1).mdtChEta()) * (tracks.at(tk2).mdtChEta()) > 0) {
1055  // check if any Ambiguity has been broken
1056  for (unsigned int AmbigTrksIdx : AmbigTrks) {
1057  if (tk2 == AmbigTrksIdx) {
1058  isResolved = true;
1059  break;
1060  }
1061  }
1062  if (isResolved) continue;
1063  // get a point on the track
1064  double Trk2ML1R = tracks.at(tk2).getML1seg().globalPosition().perp();
1065  double Trk2ML1Z = tracks.at(tk2).getML1seg().globalPosition().z();
1066  double Trk2ML2R = tracks.at(tk2).getML2seg().globalPosition().perp();
1067  double Trk2ML2Z = tracks.at(tk2).getML2seg().globalPosition().z();
1068 
1069  // find the distance between the tracks
1070  double DistML1(1000), DistML2(1000);
1071  if (tk1_isBarrel) {
1072  DistML1 = std::abs(Trk1ML1Z - Trk2ML1Z);
1073  DistML2 = std::abs(Trk1ML2Z - Trk2ML2Z);
1074  } else if (tk1_isEndcap) {
1075  DistML1 = std::abs(Trk1ML1R - Trk2ML1R);
1076  DistML2 = std::abs(Trk1ML2R - Trk2ML2R);
1077  }
1078  if (DistML1 < 40 || DistML2 < 40) {
1079  // find how many MDTs the tracks share
1080  std::vector<const Muon::MdtPrepData*> mdts1 = tracks.at(tk1).mdtHitsOnTrack();
1081  std::vector<const Muon::MdtPrepData*> mdts2 = tracks.at(tk2).mdtHitsOnTrack();
1082  nShared = 0;
1083  for (const Muon::MdtPrepData *mdt1 : mdts1) {
1084  for (const Muon::MdtPrepData *mdt2 : mdts2) {
1085  if (mdt1->identify() == mdt2->identify()) {
1086  ++nShared;
1087  break;
1088  }
1089  }
1090  }
1091 
1092  if (nShared <= 1) continue; // if the tracks share only 1 hits move to next track
1093  // store the track as ambiguous
1094  AmbigTracks.push_back(tracks.at(tk2));
1095  AmbigTrks.push_back(tk2);
1096  }
1097  } // end chamber match
1098  } // end tk2 loop
1099 
1100  if (AmbigTracks.size() == 1) {
1101  UniqueTracks.push_back(tracks.at(tk1));
1102  continue;
1103  }
1104  // Deal with any ambiguities
1105  // Barrel tracks
1106  if (tk1_isBarrel) {
1107  bool hasMomentum = tracks.at(tk1).charge() != 0;
1108  double aveX(0), aveY(0), aveZ(0), aveAlpha(0);
1109  double aveP(0), nAmbigP(0), TrkCharge(tracks.at(tk1).charge());
1110  bool allSameSign(true);
1111 
1112  for (const Tracklet &AmbigTrack : AmbigTracks) {
1113  if (!hasMomentum) {
1114  aveX += AmbigTrack.globalPosition().x();
1115  aveY += AmbigTrack.globalPosition().y();
1116  aveZ += AmbigTrack.globalPosition().z();
1117  aveAlpha += AmbigTrack.getML1seg().alpha();
1118  } else {
1119  // check the charge is the same
1120  if (std::abs(AmbigTrack.charge() - TrkCharge) > 0.1) allSameSign = false;
1121  // find the average momentum
1122  aveP += AmbigTrack.momentum().mag();
1123  ++nAmbigP;
1124  aveAlpha += AmbigTrack.alpha();
1125  aveX += AmbigTrack.globalPosition().x();
1126  aveY += AmbigTrack.globalPosition().y();
1127  aveZ += AmbigTrack.globalPosition().z();
1128  }
1129  } // end loop on ambiguous tracks
1130  if (!hasMomentum) {
1131  aveX = aveX / (double)AmbigTracks.size();
1132  aveY = aveY / (double)AmbigTracks.size();
1133  aveZ = aveZ / (double)AmbigTracks.size();
1134  Amg::Vector3D gpos(aveX, aveY, aveZ);
1135  aveAlpha = aveAlpha / (double)AmbigTracks.size();
1136  double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1137  double rErr = tracks.at(tk1).getML1seg().rError();
1138  double zErr = tracks.at(tk1).getML1seg().zError();
1139 
1140  TrackletSegment aveSegML1(m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1141  double pT = m_maxpTot * std::sin(aveSegML1.alpha());
1142  double pz = m_maxpTot * std::cos(aveSegML1.alpha());
1143  Amg::Vector3D momentum(pT * std::cos(aveSegML1.globalPosition().phi()),
1144  pT * std::sin(aveSegML1.globalPosition().phi()),
1145  pz);
1146  AmgSymMatrix(5) matrix;
1147  matrix.setIdentity();
1148  matrix(0, 0) = std::pow(tracks.at(tk1).getML1seg().rError(),2); // delta R
1149  matrix(1, 1) = std::pow(tracks.at(tk1).getML1seg().zError(),2); // delta z
1150  matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1151  matrix(3, 3) = std::pow(tracks.at(tk1).getML1seg().alphaError(),2); // delta theta
1152  matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p
1153  Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1154  UniqueTracks.push_back(aveTrack);
1155  } else if (allSameSign) {
1156  aveP = aveP / nAmbigP;
1157  double pT = aveP * std::sin(tracks.at(tk1).getML1seg().alpha());
1158  double pz = aveP * std::cos(tracks.at(tk1).getML1seg().alpha());
1159  Amg::Vector3D momentum(pT * std::cos(tracks.at(tk1).globalPosition().phi()),
1160  pT * std::sin(tracks.at(tk1).globalPosition().phi()),
1161  pz);
1162  Tracklet MyTrack = tracks.at(tk1);
1163  MyTrack.momentum(momentum);
1164  MyTrack.charge(tracks.at(tk1).charge());
1165  UniqueTracks.push_back(MyTrack);
1166  } else {
1167  aveX = aveX / (double)AmbigTracks.size();
1168  aveY = aveY / (double)AmbigTracks.size();
1169  aveZ = aveZ / (double)AmbigTracks.size();
1170  Amg::Vector3D gpos(aveX, aveY, aveZ);
1171  aveAlpha = aveAlpha / (double)AmbigTracks.size();
1172  double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1173  double rErr = tracks.at(tk1).getML1seg().rError();
1174  double zErr = tracks.at(tk1).getML1seg().zError();
1175 
1176  TrackletSegment aveSegML1(m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1177  double pT = m_maxpTot * std::sin(aveSegML1.alpha());
1178  double pz = m_maxpTot * std::cos(aveSegML1.alpha());
1179  Amg::Vector3D momentum(pT * std::cos(aveSegML1.globalPosition().phi()),
1180  pT * std::sin(aveSegML1.globalPosition().phi()),
1181  pz);
1182  AmgSymMatrix(5) matrix;
1183  matrix.setIdentity();
1184  matrix(0, 0) = std::pow(tracks.at(tk1).getML1seg().rError(),2); // delta R
1185  matrix(1, 1) = std::pow(tracks.at(tk1).getML1seg().zError(),2); // delta z
1186  matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1187  matrix(3, 3) = std::pow(tracks.at(tk1).getML1seg().alphaError(),2); // delta theta
1188  matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p
1189  Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1190  UniqueTracks.push_back(aveTrack);
1191  }
1192  } // end barrel Tracks
1193 
1194  // Endcap tracks
1195  else if (tk1_isEndcap) {
1196  std::vector<const Muon::MdtPrepData*> AllMdts;
1197  for (Tracklet const &AmbigTrack : AmbigTracks) {
1198  std::vector<const Muon::MdtPrepData*> mdts = AmbigTrack.mdtHitsOnTrack();
1199  std::vector<const Muon::MdtPrepData*> tmpAllMdt = AllMdts;
1200  for (const Muon::MdtPrepData *mdt : mdts) {
1201  bool isNewHit = true;
1202  for (const Muon::MdtPrepData *tmpmdt : tmpAllMdt) {
1203  if (mdt->identify() == tmpmdt->identify()) {
1204  isNewHit = false;
1205  break;
1206  }
1207  }
1208  if (isNewHit) AllMdts.push_back(mdt);
1209  } // end loop on mdts
1210  } // end loop on ambiguous tracks
1211 
1212  std::vector<TrackletSegment> MyECsegs = TrackletSegmentFitter(AllMdts);
1213  if (!MyECsegs.empty()) {
1214  TrackletSegment ECseg = MyECsegs.at(0);
1215  ECseg.clearMdt();
1216  double pT = m_maxpTot * std::sin(ECseg.alpha());
1217  double pz = m_maxpTot * std::cos(ECseg.alpha());
1218  Amg::Vector3D momentum(pT * std::cos(ECseg.globalPosition().phi()),
1219  pT * std::sin(ECseg.globalPosition().phi()),
1220  pz);
1221  AmgSymMatrix(5) matrix;
1222  matrix.setIdentity();
1223  matrix(0, 0) = std::pow(ECseg.rError(),2); // delta R
1224  matrix(1, 1) = std::pow(ECseg.zError(),2); // delta z
1225  matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1226  matrix(3, 3) = std::pow(ECseg.alphaError(),2); // delta theta
1227  matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
1228  Tracklet MyCombTrack(MyECsegs.at(0), ECseg, momentum, matrix, 0);
1229  UniqueTracks.push_back(MyCombTrack);
1230  } else
1231  UniqueTracks.push_back(tracks.at(tk1));
1232  } // end endcap tracks
1233 
1234  } // end loop on tracks -- tk1
1235 
1236  return UniqueTracks;
1237  }
1238 
1239 } // namespace Muon
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
base
std::string base
Definition: hcg.cxx:78
TrackletSegment::clearMdt
void clearMdt()
Definition: TrackletSegment.cxx:26
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
CXXUTILS_TRAPPING_FP
#define CXXUTILS_TRAPPING_FP
Definition: trapping_fp.h:24
Muon::MSVertexTrackletTool::m_EndcapDeltaAlphaCut
Gaudi::Property< double > m_EndcapDeltaAlphaCut
Definition: MSVertexTrackletTool.h:60
CalculateHighPtTerm.pT
pT
Definition: ICHEP2016/CalculateHighPtTerm.py:57
xAOD::identify
const Identifier & identify(const UncalibratedMeasurement *meas)
Returns the associated identifier from the muon measurement.
Definition: MuonSpectrometer/MuonPhaseII/Event/xAOD/xAODMuonPrepData/Root/UtilFunctions.cxx:61
Muon::MSVertexTrackletTool::TrackMomentum
double TrackMomentum(const Identifier trkID, const double deltaAlpha) const
Definition: MSVertexTrackletTool.cxx:920
fitman.sy
sy
Definition: fitman.py:524
Muon::MSVertexTrackletTool::SeedResiduals
static double SeedResiduals(const std::vector< const Muon::MdtPrepData * > &mdts, double slope, double inter)
Definition: MSVertexTrackletTool.cxx:569
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
add-xsec-uncert-quadrature-N.alpha
alpha
Definition: add-xsec-uncert-quadrature-N.py:110
fitman.sz
sz
Definition: fitman.py:527
Muon::MSVertexTrackletTool::m_TPContainer
SG::WriteHandleKey< xAOD::TrackParticleContainer > m_TPContainer
Definition: MSVertexTrackletTool.h:46
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
Amg::compress
void compress(const AmgSymMatrix(N) &covMatrix, std::vector< float > &vec)
Definition: EventPrimitivesHelpers.h:56
Muon::c_BMS
constexpr double c_BMS
Definition: MSVertexTrackletTool.cxx:22
Muon::MSVertexTrackletTool::CleanSegments
std::vector< TrackletSegment > CleanSegments(const std::vector< TrackletSegment > &segs) const
Definition: MSVertexTrackletTool.cxx:738
Muon::c_BOL
constexpr double c_BOL
Definition: MSVertexTrackletTool.cxx:24
Muon::MSVertexTrackletTool::m_errorCutOff
Gaudi::Property< double > m_errorCutOff
Definition: MSVertexTrackletTool.h:54
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
Muon::MSVertexTrackletTool::m_straightTrackletInvPerr
Gaudi::Property< double > m_straightTrackletInvPerr
Definition: MSVertexTrackletTool.h:65
hist_file_dump.d
d
Definition: hist_file_dump.py:137
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
skel.it
it
Definition: skel.GENtoEVGEN.py:396
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Muon::c_BML
constexpr double c_BML
Definition: MSVertexTrackletTool.cxx:23
Muon::MSVertexTrackletTool::m_straightTrackletpTot
Gaudi::Property< double > m_straightTrackletpTot
Definition: MSVertexTrackletTool.h:64
xAOD::TrackParticle_v1::setDefiningParameters
void setDefiningParameters(float d0, float z0, float phi0, float theta, float qOverP)
Set the defining parameters.
Definition: TrackParticle_v1.cxx:177
Trk::z0
@ z0
Definition: ParamDefs.h:64
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
TrackletSegment::getChMidPoint
double getChMidPoint() const
Definition: TrackletSegment.cxx:44
Trk::locR
@ locR
Definition: ParamDefs.h:44
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Muon::MSVertexTrackletTool::initialize
virtual StatusCode initialize() override
Definition: MSVertexTrackletTool.cxx:34
Muon::MSVertexTrackletTool::findTracklets
StatusCode findTracklets(std::vector< Tracklet > &tracklets, const EventContext &ctx) const override
Definition: MSVertexTrackletTool.cxx:44
Muon
NRpcCablingAlg reads raw condition data and writes derived condition data to the condition store.
Definition: TrackSystemController.h:45
Muon::MSVertexTrackletTool::MSVertexTrackletTool
MSVertexTrackletTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: MSVertexTrackletTool.cxx:27
empty
bool empty(TH1 *h)
Definition: computils.cxx:295
Muon::MSVertexTrackletTool::DeltabCalc
bool DeltabCalc(const TrackletSegment &ML1seg, const TrackletSegment &ML2seg) const
Definition: MSVertexTrackletTool.cxx:887
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
Muon::MdtStatusDriftTime
@ MdtStatusDriftTime
The tube produced a vaild measurement.
Definition: MdtDriftCircleStatus.h:34
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
TrackletSegment::getIdentifier
const Identifier getIdentifier() const
Definition: TrackletSegment.cxx:31
Muon::MSVertexTrackletTool::convertToTrackParticles
static void convertToTrackParticles(std::vector< Tracklet > &tracklets, SG::WriteHandle< xAOD::TrackParticleContainer > &container)
Definition: MSVertexTrackletTool.cxx:330
Muon::MSVertexTrackletTool::TrackletSegmentFitterCore
std::vector< TrackletSegment > TrackletSegmentFitterCore(const std::vector< const Muon::MdtPrepData * > &mdts, const std::vector< std::pair< double, double > > &SeedParams) const
Definition: MSVertexTrackletTool.cxx:584
TrackletSegment::rError
double rError() const
Definition: TrackletSegment.cxx:36
Muon::MSVertexTrackletTool::m_d12_max
Gaudi::Property< double > m_d12_max
Definition: MSVertexTrackletTool.h:52
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition: Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
TrackletSegment::globalPosition
const Amg::Vector3D & globalPosition() const
Definition: TrackletSegment.cxx:32
xAOD::TrackParticle_v1::setFitQuality
void setFitQuality(float chiSquared, float numberDoF)
Set the 'Fit Quality' information.
Definition: TrackParticle_v1.cxx:534
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
DiTauMassTools::ignore
void ignore(T &&)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:58
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
id2
HWIdentifier id2
Definition: LArRodBlockPhysicsV0.cxx:562
lumiFormat.i
int i
Definition: lumiFormat.py:85
TrackletSegment
New segment class for single ML segments.
Definition: TrackletSegment.h:20
Muon::MSVertexTrackletTool::TrackMomentumError
double TrackMomentumError(const TrackletSegment &ml1, const TrackletSegment &ml2) const
Definition: MSVertexTrackletTool.cxx:950
TrackletSegment::alphaError
double alphaError() const
Definition: TrackletSegment.cxx:34
Trk::theta
@ theta
Definition: ParamDefs.h:66
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
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
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
TrackletSegment::alpha
double alpha() const
Definition: TrackletSegment.cxx:33
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:523
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
xAOD::LowPtTrack
@ LowPtTrack
A LowPt track.
Definition: TrackingPrimitives.h:77
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Muon::MSVertexTrackletTool::m_SeedResidual
Gaudi::Property< double > m_SeedResidual
Definition: MSVertexTrackletTool.h:56
Muon::MSVertexTrackletTool::m_minSegFinderChi2
Gaudi::Property< double > m_minSegFinderChi2
Definition: MSVertexTrackletTool.h:57
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
Muon::MSVertexTrackletTool::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MSVertexTrackletTool.h:48
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
Tracklet::charge
void charge(double charge)
Definition: Tracklet.cxx:41
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
Muon::MSVertexTrackletTool::IgnoreMDTChamber
bool IgnoreMDTChamber(const Muon::MdtPrepData *mdtHit) const
Definition: MSVertexTrackletTool.cxx:356
Muon::MuonPrepDataCollection
Template to hold collections of MuonPrepRawData objects.
Definition: MuonPrepDataCollection.h:46
xAOD::TrackParticle_v1::setTrackProperties
void setTrackProperties(const TrackProperties properties)
Methods setting the TrackProperties.
Muon::MSVertexTrackletTool::m_mdtTESKey
SG::ReadHandleKey< Muon::MdtPrepDataContainer > m_mdtTESKey
Definition: MSVertexTrackletTool.h:45
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
Muon::MSVertexTrackletTool::m_maxpTot
Gaudi::Property< double > m_maxpTot
Definition: MSVertexTrackletTool.h:63
trapping_fp.h
Tell the compiler to optimize assuming that FP may trap.
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
Trk::d0
@ d0
Definition: ParamDefs.h:63
Muon::MSVertexTrackletTool::SortMDThits
int SortMDThits(std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt, const EventContext &ctx) const
Definition: MSVertexTrackletTool.cxx:380
Amg::error
double error(const Amg::MatrixX &mat, int index)
return diagonal error of the matrix caller should ensure the matrix is symmetric and the index is in ...
Definition: EventPrimitivesHelpers.h:40
charge
double charge(const T &p)
Definition: AtlasPID.h:756
Muon::MSVertexTrackletTool::addMDTHits
void addMDTHits(std::vector< const Muon::MdtPrepData * > &hits, std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt) const
Definition: MSVertexTrackletTool.cxx:432
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
Tracklet::momentum
void momentum(const Amg::Vector3D &p)
Definition: Tracklet.cxx:40
Muon::MSVertexTrackletTool::m_d13_max
Gaudi::Property< double > m_d13_max
Definition: MSVertexTrackletTool.h:53
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
library_scraper.dd
list dd
Definition: library_scraper.py:46
Muon::MSVertexTrackletTool::m_BarrelDeltaAlphaCut
Gaudi::Property< double > m_BarrelDeltaAlphaCut
Definition: MSVertexTrackletTool.h:59
Muon::c_BIL
constexpr double c_BIL
Definition: MSVertexTrackletTool.cxx:21
xAOD::TrackParticle_v1::setDefiningParametersCovMatrixVec
void setDefiningParametersCovMatrixVec(const std::vector< float > &cov)
Definition: TrackParticle_v1.cxx:460
TrackletSegment::zError
double zError() const
Definition: TrackletSegment.cxx:35
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:76
python.testIfMatch.matrix
matrix
Definition: testIfMatch.py:66
MSVertexTrackletTool.h
python.compareTCTs.ratio
ratio
Definition: compareTCTs.py:295
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
Muon::MSVertexTrackletTool::m_maxDeltabCut
Gaudi::Property< double > m_maxDeltabCut
Definition: MSVertexTrackletTool.h:61
DEBUG
#define DEBUG
Definition: page_access.h:11
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
Muon::MSVertexTrackletTool::TrackletSegmentFitter
std::vector< TrackletSegment > TrackletSegmentFitter(const std::vector< const Muon::MdtPrepData * > &mdts) const
Definition: MSVertexTrackletTool.cxx:453
python.LArCondContChannels.isBarrel
isBarrel
Definition: LArCondContChannels.py:659
Muon::MSVertexTrackletTool::m_minpTot
Gaudi::Property< double > m_minpTot
Definition: MSVertexTrackletTool.h:62
Muon::MSVertexTrackletTool::ResolveAmbiguousTracklets
std::vector< Tracklet > ResolveAmbiguousTracklets(std::vector< Tracklet > &tracks) const
Definition: MSVertexTrackletTool.cxx:1008
Muon::MSVertexTrackletTool::m_tightTrackletRequirement
Gaudi::Property< bool > m_tightTrackletRequirement
Definition: MSVertexTrackletTool.h:66
xAOD::track
@ track
Definition: TrackingPrimitives.h:512
Muon::MSVertexTrackletTool::SegSeeds
std::vector< std::pair< double, double > > SegSeeds(const std::vector< const Muon::MdtPrepData * > &mdts) const
Definition: MSVertexTrackletTool.cxx:465
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Tracklet
Definition: Tracklet.h:17
Trk::phi0
@ phi0
Definition: ParamDefs.h:65
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Fri Jan 10 2025 21:14:29 for ATLAS Offline Software by   doxygen 1.8.18