80 ATH_CHECK(container.record(std::make_unique<xAOD::TrackParticleContainer>(), std::make_unique<xAOD::TrackParticleAuxContainer>()));
83 std::vector<std::vector<const Muon::MdtPrepData*> > SortedMdt;
87 if (nMDT <= 0) {
return StatusCode::SUCCESS; }
100 double d12_max = 50.;
101 double d13_max = 80.;
103 constexpr
double errorCutOff = 0.001;
104 std::vector<TrackletSegment> segs[6][2][16];
105 std::vector<std::vector<const Muon::MdtPrepData*> >::const_iterator ChamberItr = SortedMdt.begin();
106 for (; ChamberItr != SortedMdt.end(); ++ChamberItr) {
107 std::vector<TrackletSegment> mlsegments;
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;
114 if (stName == 1 && std::abs(stEta) >= 7)
continue;
115 if (stName == 53 || stName == 54)
continue;
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;
121 int maxLayer =
m_idHelperSvc->mdtIdHelper().tubeLayerMax((*mdt1)->identify());
122 int ML =
m_idHelperSvc->mdtIdHelper().multilayer((*mdt1)->identify());
123 for (; mdt1 != mdtEnd; ++mdt1) {
129 int tl1 =
m_idHelperSvc->mdtIdHelper().tubeLayer((*mdt1)->identify());
130 if (tl1 == maxLayer)
break;
131 std::vector<const Muon::MdtPrepData*>::const_iterator mdt2 = (mdt1 + 1);
132 for (; mdt2 != mdtEnd; ++mdt2) {
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))
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;
150 if ((tl2 - tl1) == 0 && (
m_idHelperSvc->mdtIdHelper().tube((*mdt2)->identify()) -
154 std::vector<const Muon::MdtPrepData*>::const_iterator mdt3 = (mdt2 + 1);
155 for (; mdt3 != mdtEnd; ++mdt3) {
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()) -
169 if (std::abs((*mdt3)->globalPosition().z() - (*mdt1)->globalPosition().z()) > d13_max &&
170 (stName <= 11 || stName == 52))
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;
179 std::vector<const Muon::MdtPrepData*> mdts;
180 mdts.push_back((*mdt1));
181 mdts.push_back((*mdt2));
182 mdts.push_back((*mdt3));
184 for (
unsigned int i = 0;
i < tmpSegs.size(); ++
i) mlsegments.push_back(tmpSegs.at(
i));
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));
206 ATH_MSG_WARNING(
"Found segments belonging to chamber " << stName <<
" that have not been stored");
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]);
221 for (
int st = 0; st < 6; ++st) {
224 for (
int sector = 0; sector < 16; ++sector) {
225 for (
unsigned int i1 = 0; i1 < CleanSegs[st][0][sector].size(); ++i1) {
227 int stName = CleanSegs[st][0][sector].at(i1).mdtChamber();
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;
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())
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));
258 if (std::abs(deltaAlpha) < DeltaAlphaCut && goodDeltab) {
261 if (deltaAlpha * CleanSegs[st][0][sector].at(i1).globalPosition().
z() *
262 std::tan(CleanSegs[st][0][sector].at(i1).alpha()) <
265 float pTot =
TrackMomentum(CleanSegs[st][0][sector].at(i1).mdtChamber(), deltaAlpha);
266 if (pTot < 800.)
continue;
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));
274 if (!CombinedSeg.empty()) {
277 float pT = pTot *
std::sin(CombinedSeg[0].alpha());
278 float pz = pTot *
std::cos(CombinedSeg[0].alpha());
284 matrix(0, 0) =
sq(CombinedSeg[0].rError());
285 matrix(1, 1) =
sq(CombinedSeg[0].zError());
288 matrix(3, 3) =
sq(CombinedSeg[0].alphaError());
294 tracklets.push_back(tmpTrk);
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());
303 pT *
std::sin(CleanSegs[st][0][sector].at(i1).globalPosition().
phi()),
pz);
307 matrix(0, 0) =
sq(CleanSegs[st][0][sector].at(i1).rError());
308 matrix(1, 1) =
sq(CleanSegs[st][0][sector].at(i1).zError());
311 matrix(3, 3) =
sq(CleanSegs[st][0][sector].at(i1).alphaError());
318 tracklets.push_back(tmpTrk);
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));
327 if (!CombinedSeg.empty()) {
329 float pT = 100000.0 *
std::sin(CombinedSeg[0].alpha());
330 float pz = 100000.0 *
std::cos(CombinedSeg[0].alpha());
334 matrix(0, 0) =
sq(CombinedSeg[0].rError());
335 matrix(1, 1) =
sq(CombinedSeg[0].zError());
338 matrix(3, 3) =
sq(CombinedSeg[0].alphaError());
342 pT * std::
sin(CombinedSeg[0].globalPosition().phi()),
pz);
344 tracklets.push_back(tmpTrk);