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