MuonTrackSummaryHelperTool.cxx 26.9 KB
Newer Older
1
/*
2
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3
4
*/

5
#include "MuonTrackSummaryHelperTool.h"
6

7
8
9
10
#include <cassert>
#include <cmath>
#include <set>

11
#include "MuonCompetingRIOsOnTrack/CompetingMuonClustersOnTrack.h"
12
#include "MuonPrepRawData/CscClusterStatus.h"
13
#include "MuonRIO_OnTrack/CscClusterOnTrack.h"
14
#include "MuonRIO_OnTrack/MdtDriftCircleOnTrack.h"
15
#include "MuonRIO_OnTrack/MMClusterOnTrack.h"
16
#include "MuonRIO_OnTrack/MuonDriftCircleErrorStrategy.h"
17
18
#include "MuonReadoutGeometry/MdtReadoutElement.h"
#include "MuonReadoutGeometry/MuonDetectorManager.h"
19
20
21
22
#include "StoreGate/ReadHandle.h"
#include "TrkCompetingRIOsOnTrack/CompetingRIOsOnTrack.h"
#include "TrkDetElementBase/TrkDetElementBase.h"
#include "TrkMeasurementBase/MeasurementBase.h"
23
#include "TrkParameters/TrackParameters.h"
24
#include "TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h"
25
#include "TrkSurfaces/Surface.h"
26
27
#include "TrkTrack/Track.h"
#include "TrkTrack/TrackStateOnSurface.h"
28

29
Muon::MuonTrackSummaryHelperTool::MuonTrackSummaryHelperTool(const std::string& t, const std::string& n, const IInterface* p) :
30
31
    base_class(t, n, p) {
    declareInterface<ITrackSummaryHelperTool>(this);
32
33
}

34
35
36
37
38
39
40
41
42
43
StatusCode Muon::MuonTrackSummaryHelperTool::initialize() {
    ATH_CHECK(m_DetectorManagerKey.initialize());
    if (m_calculateCloseHits && !m_extrapolator.empty()) {
        ATH_CHECK(m_extrapolator.retrieve());
    } else {
        m_extrapolator.disable();
    }
    ATH_CHECK(m_idHelperSvc.retrieve());
    ATH_CHECK(m_mdtKey.initialize());
    return StatusCode::SUCCESS;
44
45
}

46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
void Muon::MuonTrackSummaryHelperTool::analyse(const Trk::Track& /**trk*/, const Trk::RIO_OnTrack* rot,
                                               const Trk::TrackStateOnSurface* tsos, std::vector<int>& information,
                                               std::bitset<Trk::numberOfDetectorTypes>& /**hitPattern*/) const {
    using namespace Trk;
    if (tsos->type(Trk::TrackStateOnSurface::Outlier)) return;  // ignore outliers

    Identifier id = rot->identify();
    ATH_MSG_DEBUG("Processing rot: " << m_idHelperSvc->toString(id));
    if (m_idHelperSvc->isRpc(id)) {
        if (m_idHelperSvc->rpcIdHelper().measuresPhi(id))
            increment(information[numberOfRpcPhiHits]);
        else
            increment(information[numberOfRpcEtaHits]);
    } else if (m_idHelperSvc->isCsc(id)) {
        if (m_idHelperSvc->cscIdHelper().measuresPhi(id))
            increment(information[numberOfCscPhiHits]);
        else {
            increment(information[numberOfCscEtaHits]);
            const CscClusterOnTrack* clus = dynamic_cast<const CscClusterOnTrack*>(rot);
            if (clus && ((clus->status() == Muon::CscStatusUnspoiled) || (clus->status() == Muon::CscStatusSplitUnspoiled)))
                increment(information[numberOfCscUnspoiltEtaHits]);
        }
    } else if (m_idHelperSvc->isTgc(id)) {
        if (m_idHelperSvc->tgcIdHelper().isStrip(id))
            increment(information[numberOfTgcPhiHits]);
        else
            increment(information[numberOfTgcEtaHits]);
    } else if (m_idHelperSvc->isMdt(id)) {
        increment(information[numberOfMdtHits]);
    } else if (m_idHelperSvc->issTgc(id)) {
        if (m_idHelperSvc->stgcIdHelper().measuresPhi(id)) increment(information[numberOfStgcPhiHits]);
        // we do not discriminate between pads or wires
        else
            increment(information[numberOfStgcEtaHits]);
    } else if (m_idHelperSvc->isMM(id)) {
        increment(information[numberOfMmHits]);
    } else {
        ATH_MSG_ERROR("Unknown muon detector type ");
        ATH_MSG_ERROR("Dumping TrackStateOnSurface " << *tsos);
85
    }
86
    return;
87
88
}

89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
void Muon::MuonTrackSummaryHelperTool::analyse(const Trk::Track& trk, const Trk::CompetingRIOsOnTrack* crot,
                                               const Trk::TrackStateOnSurface* tsos, std::vector<int>& information,
                                               std::bitset<Trk::numberOfDetectorTypes>& hitPattern) const {
    // For competing ROTs we *only* count hits that are on different layers.
    std::set<Identifier> layIds;
    for (unsigned int i = 0; i < crot->numberOfContainedROTs(); i++) {
        const Trk::RIO_OnTrack* rot = &crot->rioOnTrack(i);
        Identifier layId = m_idHelperSvc->layerId(rot->identify());
        ATH_MSG_DEBUG("ROT " << i << "\t LayerId=" << m_idHelperSvc->toString(layId));
        std::pair<std::set<Identifier>::iterator, bool> pr = layIds.insert(layId);
        if (pr.second == true) {
            // layer not seen before
            ATH_MSG_DEBUG("Have found hit on new layer. # of layers for this cROT currently=" << layIds.size());
            analyse(trk, rot, tsos, information, hitPattern);
        }
104
105
106
    }
}

107
108
109
110
111
112
void Muon::MuonTrackSummaryHelperTool::increment(int& type) const {
    if (type < 0)
        type = 1;  // they all start off at -1, so can't just increment
    else
        ++type;
    return;
113
114
}

115
116
117
void Muon::MuonTrackSummaryHelperTool::searchForHoles(const Trk::Track& /**track*/, std::vector<int>& /**information*/,
                                                      Trk::ParticleHypothesis /**hyp*/) const {
    ATH_MSG_WARNING("searchForHoles is not implemented in MuonTrackSummaryHelperTool");
118
119
}

120
121
122
123
124
void Muon::MuonTrackSummaryHelperTool::addDetailedTrackSummary(const Trk::Track& track, Trk::TrackSummary& summary) const {
    if (summary.m_muonTrackSummary) {
        ATH_MSG_DEBUG("TrackSummary already has detailed muon track summary, not adding a new one");
        return;
    }
125

126
127
128
129
130
131
    SG::ReadCondHandle<MuonGM::MuonDetectorManager> DetectorManagerHandle{m_DetectorManagerKey};
    const MuonGM::MuonDetectorManager* MuonDetMgr{*DetectorManagerHandle};
    if (MuonDetMgr == nullptr) {
        ATH_MSG_ERROR("Null pointer to the read MuonDetectorManager conditions object");
        return;
    }
132

133
134
135
136
137
    ATH_MSG_DEBUG("Adding detailed muon track summary");
    ATH_MSG_DEBUG(track.info());
    // loop over track and get chamber Identifiers
    const DataVector<const Trk::TrackStateOnSurface>* states = track.trackStateOnSurfaces();
    if (!states || states->empty()) { return; }
138

139
140
    Trk::MuonTrackSummary* muonTrackSummary = new Trk::MuonTrackSummary();
    Trk::MuonTrackSummary& trackSummary = *muonTrackSummary;
141

142
143
    Trk::MuonTrackSummary::ChamberHitSummary* currentChamberSummary = nullptr;
    const Trk::TrackParameters* currentChamberPars = nullptr;
144

145
146
147
148
149
    // loop over TSOSs
    DataVector<const Trk::TrackStateOnSurface>::const_iterator tsit = states->begin();
    DataVector<const Trk::TrackStateOnSurface>::const_iterator tsit_end = states->end();
    for (; tsit != tsit_end; ++tsit) {
        const Trk::TrackParameters* pars = (*tsit)->trackParameters();
150

151
152
153
154
        if ((*tsit)->type(Trk::TrackStateOnSurface::Scatterer)) {
            ++trackSummary.m_nscatterers;
            continue;
        }
155

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
        if ((*tsit)->type(Trk::TrackStateOnSurface::Hole)) {
            if (!pars) {
                ATH_MSG_WARNING(" Hole state without track parameters, cannot identify hole ");
                continue;
            }
            if (pars->associatedSurface().associatedDetectorElement()) {
                Identifier id = pars->associatedSurface().associatedDetectorElement()->identify();
                bool issTgc = m_idHelperSvc->issTgc(id);
                if (issTgc) {
                    // get the identifier for phi or eta holes
                    Identifier idh = pars->associatedSurface().associatedDetectorElementIdentifier();
                    if (idh.is_valid()) { id = idh; }
                }
                if (!id.is_valid() || !m_idHelperSvc->isMuon(id)) continue;
                Identifier chId = m_idHelperSvc->chamberId(id);
                // for is summary sTGC split STGC1 and STGC2
                if (issTgc) chId = m_idHelperSvc->detElId(id);
                bool isFirst = isFirstProjection(id);
                bool isMdt = m_idHelperSvc->isMdt(id);

                // check whether first chamber or new chamber
                if (!currentChamberSummary || currentChamberSummary->m_chId != chId) {
                    /** go back through the chamber, nallhitsinroad - nhitsontrack = ncloseHits
                        for both projection1 and projection2
                        pass chamber summary to function calculateRoadHits(currentChamberSummary) */
                    if (m_calculateCloseHits && currentChamberSummary && currentChamberPars) {
                        ATH_MSG_VERBOSE(" Calculating close hits (last hit a hole)");
                        calculateRoadHits(*currentChamberSummary, *currentChamberPars);
                    }

                    // given that we cannot separate eta/phi holes, redo the assignment before moving to the next chamber
                    if (currentChamberSummary && !currentChamberSummary->isMdt()) { updateHoleContent(*currentChamberSummary); }

                    ATH_MSG_VERBOSE(" Adding new chamber (holes) " << m_idHelperSvc->toString(id) << " " << *pars);
                    trackSummary.m_chamberHitSummary.push_back(Trk::MuonTrackSummary::ChamberHitSummary(chId, isMdt));
                    currentChamberSummary = &trackSummary.m_chamberHitSummary.back();
                    currentChamberPars = pars;
                }
                if (!issTgc) {
                    Trk::MuonTrackSummary::ChamberHitSummary::Projection& proj =
                        isFirst ? currentChamberSummary->m_first : currentChamberSummary->m_second;
                    ++proj.nholes;
                } else {
                    // sTgc holes keep track of phi and eta
                    if (m_idHelperSvc->measuresPhi(id)) {
                        ATH_MSG_VERBOSE(" counting sTGC phi hole ");
                        Trk::MuonTrackSummary::ChamberHitSummary::Projection& proj = currentChamberSummary->m_second;
                        ++proj.nholes;
                    } else {
                        ATH_MSG_VERBOSE(" counting sTGC eta hole ");
                        Trk::MuonTrackSummary::ChamberHitSummary::Projection& proj = currentChamberSummary->m_first;
                        ++proj.nholes;
                    }
                }
            }
            continue;
        }
213

214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
        // check whether state is a measurement
        const Trk::MeasurementBase* meas = (*tsit)->measurementOnTrack();
        if (!meas) { continue; }

        const Trk::PseudoMeasurementOnTrack* pseudo = dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(meas);
        if (pseudo) {
            ++trackSummary.m_npseudoMeasurements;
            continue;
        }

        if (!pars) {
            ATH_MSG_DEBUG("measurement without pars");
            continue;
        }

        Amg::Vector2D locPos;
        if (!meas->associatedSurface().globalToLocal(pars->position(), pars->position(), locPos)) {
            ATH_MSG_DEBUG(" localToGlobal failed !!!!! ");
            continue;
        }
        bool inBounds = true;

        Identifier id;
        std::set<Identifier> layIds;
        std::set<Identifier> goodLayIds;  // holds mdt hits that have not been deweighted

        // check whether ROT
        const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(meas);
        if (rot) {
            id = rot->identify();
            if (!m_idHelperSvc->isMuon(id)) continue;

            // bound checks
            double tol1 = 100.;
            double tol2 = 2 * tol1;
            if (!pseudo && m_idHelperSvc->isMdt(id)) tol1 = 5.;
250

251
252
            // we need a special bound check for MDTs so we cast to SL surface
            const Trk::StraightLineSurface* slSurf = dynamic_cast<const Trk::StraightLineSurface*>(&meas->associatedSurface());
253
254
255
            // we need a special bound check also for MMs to consider edge passivation
            const MMClusterOnTrack* mmClusterOnTrack = dynamic_cast<const MMClusterOnTrack*>(meas);

256
257
258
            if (slSurf) {
                // perform bound check only for second coordinate
                inBounds = slSurf->bounds().insideLoc2(locPos, tol2);
259
260
261
            } else if (mmClusterOnTrack) {
                // for MM, perform the bound check from the detector element
                inBounds = mmClusterOnTrack->detectorElement()->insideActiveBounds(id, locPos, tol1, tol2);
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
            } else {
                inBounds = meas->associatedSurface().insideBounds(locPos, tol1, tol2);
            }

            Identifier layId = m_idHelperSvc->layerId(id);
            layIds.insert(layId);
            const MdtDriftCircleOnTrack* mdtdc = dynamic_cast<const MdtDriftCircleOnTrack*>(rot);
            if (mdtdc) {
                MuonDriftCircleErrorStrategy errStrat = mdtdc->errorStrategy();
                if (!errStrat.creationParameter(MuonDriftCircleErrorStrategy::FixedError) &&
                    !errStrat.creationParameter(MuonDriftCircleErrorStrategy::StationError)) {
                    goodLayIds.insert(layId);
                }
            } else if (m_idHelperSvc->isCsc(id)) {
                const Muon::CscClusterOnTrack* cscClus = dynamic_cast<const Muon::CscClusterOnTrack*>(rot);
                if (cscClus->status() == Muon::CscClusterStatus::CscStatusUnspoiled ||
                    cscClus->status() == Muon::CscClusterStatus::CscStatusSplitUnspoiled)
                    goodLayIds.insert(layId);
            } else if (m_idHelperSvc->isMM(id)) {
                // MM quality requirements to be inserted here if needed
                goodLayIds.insert(layId);
            } else if (m_idHelperSvc->issTgc(id)) {
                // sTGC quality requirements to be inserted here if needed
                goodLayIds.insert(layId);
            }
        } else {
            const Muon::CompetingMuonClustersOnTrack* crot = dynamic_cast<const Muon::CompetingMuonClustersOnTrack*>(meas);
            if (crot) {
                if (crot->containedROTs().empty()) continue;

                // take id of first ROT
                id = crot->containedROTs().front()->identify();

                // count layers in competing rot
                std::vector<const Muon::MuonClusterOnTrack*>::const_iterator cl_it = crot->containedROTs().begin();
                std::vector<const Muon::MuonClusterOnTrack*>::const_iterator cl_it_end = crot->containedROTs().end();
                for (; cl_it != cl_it_end; ++cl_it) {
                    // get layer Identifier and insert it into set
                    Identifier layId = m_idHelperSvc->layerId((*cl_it)->identify());
                    layIds.insert(layId);
                    if (m_idHelperSvc->isCsc(id)) {
                        const Muon::CscClusterOnTrack* cscClus = dynamic_cast<const Muon::CscClusterOnTrack*>(rot);
                        if (cscClus->status() == Muon::CscClusterStatus::CscStatusUnspoiled ||
                            cscClus->status() == Muon::CscClusterStatus::CscStatusSplitUnspoiled)
                            goodLayIds.insert(layId);
                    }
                }
            } else {
                continue;
            }
312
        }
313
        Identifier chId = m_idHelperSvc->chamberId(id);
314
315
316
317
        // for is summary sTGC split STGC1 and STGC2
        bool issTgc = m_idHelperSvc->issTgc(id);
        if (issTgc) chId = m_idHelperSvc->detElId(id);
        bool isFirst = isFirstProjection(id);
318
        bool isMdt = m_idHelperSvc->isMdt(id);
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
        ATH_MSG_VERBOSE(" Adding hit " << m_idHelperSvc->toString(id));

        /** check whether first chamber or new chamber */
        if (!currentChamberSummary || currentChamberSummary->m_chId != chId) {
            if (m_calculateCloseHits && currentChamberSummary && currentChamberPars) {
                ATH_MSG_VERBOSE(" Calculating close hits (last hit a measurement)");
                calculateRoadHits(*currentChamberSummary, *currentChamberPars);
            }

            // given that we cannot separate eta/phi holes, redo the assignment before moving to the next chamber
            if (currentChamberSummary && !currentChamberSummary->isMdt()) { updateHoleContent(*currentChamberSummary); }

            ATH_MSG_VERBOSE(" Adding new chamber " << m_idHelperSvc->toString(id) << " " << *pars);
            trackSummary.m_chamberHitSummary.push_back(Trk::MuonTrackSummary::ChamberHitSummary(chId, isMdt));
            currentChamberSummary = &trackSummary.m_chamberHitSummary.back();
            currentChamberPars = pars;
        }

        Trk::MuonTrackSummary::ChamberHitSummary::Projection& proj =
            isFirst ? currentChamberSummary->m_first : currentChamberSummary->m_second;

        if ((*tsit)->type(Trk::TrackStateOnSurface::Outlier)) {
            // MDTs: count outlier as delta electron if rDrift < rTrack < innerTubeRadius
            if (isMdt && pars) {
                double rDrift = std::abs(meas->localParameters()[Trk::locR]);
                double rTrack = std::abs(pars->parameters()[Trk::locR]);
                double innerRadius = MuonDetMgr->getMdtReadoutElement(id)->innerTubeRadius();
                if (rTrack > rDrift && rTrack < innerRadius) {
                    ++proj.ndeltas;
                    continue;
                }
            }
            ++proj.noutliers;
352

353
        } else {
354
355
            proj.nhits += layIds.size();
            proj.ngoodHits += goodLayIds.size();
356
        }
357
        if (!inBounds && isMdt) proj.noutBounds++;
358

359
    }  // end of for loop over Track State on Surfaces
360

361
362
363
364
365
    /** calculate road hits for last chamber on track
    (otherwise it would have been skipped) */
    if (m_calculateCloseHits && currentChamberSummary && currentChamberPars) {
        ATH_MSG_VERBOSE(" Calculating close hits (end of hit list)");
        calculateRoadHits(*currentChamberSummary, *currentChamberPars);
366
367
    }

368
369
370
    // given that we cannot separate eta/phi holes, redo the assignment before moving to the next chamber
    if (currentChamberSummary && !currentChamberSummary->isMdt()) { updateHoleContent(*currentChamberSummary); }

371
    summary.m_muonTrackSummary.reset(muonTrackSummary);
372
373
374
375
376
}

void Muon::MuonTrackSummaryHelperTool::updateHoleContent(Trk::MuonTrackSummary::ChamberHitSummary& chamberHitSummary) const {
    if (m_idHelperSvc->issTgc(chamberHitSummary.chamberId())) {
        ATH_MSG_DEBUG(" holes eta " << chamberHitSummary.etaProjection().nholes << " phi " << chamberHitSummary.phiProjection().nholes);
377
378
    }

379
380
381
382
383
384
385
    if (m_idHelperSvc->issTgc(chamberHitSummary.chamberId()) || m_idHelperSvc->isMM(chamberHitSummary.chamberId())) { return; }

    SG::ReadCondHandle<MuonGM::MuonDetectorManager> DetectorManagerHandle{m_DetectorManagerKey};
    const MuonGM::MuonDetectorManager* MuonDetMgr{*DetectorManagerHandle};
    if (!MuonDetMgr) {
        ATH_MSG_ERROR("Null pointer to the read MuonDetectorManager conditions object");
        return;
386
387
    }

388
389
390
391
392
393
394
395
396
397
398
399
    bool isCsc = m_idHelperSvc->isCsc(chamberHitSummary.chamberId());
    int neta = isCsc ? 4 : 2;
    int nphi = isCsc ? 4 : 2;
    if (m_idHelperSvc->isTgc(chamberHitSummary.chamberId())) {
        const MuonGM::TgcReadoutElement* detEl = MuonDetMgr->getTgcReadoutElement(chamberHitSummary.chamberId());
        if (!detEl) {
            ATH_MSG_WARNING(" No detector element found for " << m_idHelperSvc->toStringChamber(chamberHitSummary.chamberId()));
            return;
        }

        // get list of layers with a hole
        neta = detEl->Ngasgaps();
400
401
    }

402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
    // code to recalculate the hole counts as they are not correct.
    // This is due to the fact that the identification of the layers goes via the readout element identifier
    // so it is not possible to separate eta and phi holes
    int nMisEta = neta - chamberHitSummary.etaProjection().nhits - chamberHitSummary.etaProjection().noutliers;
    int nMisPhi = nphi - chamberHitSummary.phiProjection().nhits - chamberHitSummary.phiProjection().noutliers;
    int nholes = chamberHitSummary.etaProjection().nholes + chamberHitSummary.phiProjection().nholes;
    if (nMisEta > 0 && nholes > 0) {
        chamberHitSummary.m_first.nholes = nMisEta;
        nholes -= nMisEta;
    }
    if (nMisPhi > 0 && nholes > 0) {
        chamberHitSummary.m_second.nholes = nholes;
        if (nholes != nMisPhi) {
            ATH_MSG_DEBUG("Inconsistent hole count: expected hits eta "
                          << neta << " phi " << nphi << " hits eta "
                          << chamberHitSummary.etaProjection().nhits + chamberHitSummary.etaProjection().noutliers << " phi "
                          << chamberHitSummary.phiProjection().nhits + chamberHitSummary.phiProjection().noutliers << " missed eta "
                          << nMisEta << " phi " << nMisPhi << " holes eta " << chamberHitSummary.etaProjection().nholes << " phi "
                          << chamberHitSummary.phiProjection().nholes << " recalculated phi holes " << nholes);
421
        }
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
    }
}

void Muon::MuonTrackSummaryHelperTool::calculateRoadHits(Trk::MuonTrackSummary::ChamberHitSummary& chamberHitSummary,
                                                         const Trk::TrackParameters& pars) const {
    bool isStraightLine = false;
    if (pars.parameters().rows() < 5) {  // no momentum parameter given
        isStraightLine = true;
    } else if (std::abs(pars.parameters()[4]) < 1e-8) {  // |p| > 1e8 MeV = 100 TeV
        isStraightLine = true;
    } else {
        // Determine if TrackParameters correspond to a straight track (the ugly way)
        if (pars.covariance()) {
            const AmgSymMatrix(5)& covMat = *pars.covariance();
            if (covMat.rows() < 5) {  // no momentum available
                isStraightLine = true;
            } else {
                // if no error on momentum given, assume no momentum was measured (extrapolator fails on zero error)
                if (std::abs(covMat(4, 4)) < 1e-20) isStraightLine = true;
            }
442
        }
443
    }
444
445
446
447
448
449
450
    const Trk::IExtrapolator* extrapolator = 0;
    if (!isStraightLine && !m_extrapolator.empty()) {
        extrapolator = &*(m_extrapolator);
    } else if (isStraightLine && !m_slExtrapolator.empty()) {
        extrapolator = &*(m_slExtrapolator);
    }
    if (!extrapolator) return;
451

452
    ATH_MSG_DEBUG("road hits for chamber " << m_idHelperSvc->toString(chamberHitSummary.chamberId()));
453

454
455
    // currently treating MDTs only
    if (!chamberHitSummary.isMdt()) return;
456

457
458
459
460
461
462
    // loop over Mdt Prds (all hits in this chamber) and try to find prds of tubes with hits
    const Muon::MdtPrepDataCollection* mdtPrdCol = findMdtPrdCollection(chamberHitSummary.chamberId());
    if (!mdtPrdCol) {
        ATH_MSG_DEBUG(" Retrieval of MdtPrepDataCollection failed!! ");
        return;
    }
463

464
465
466
467
468
    if (isStraightLine) {
        ATH_MSG_VERBOSE("Doing straight line extrapolation to get hits in road");
    } else {
        ATH_MSG_VERBOSE("Doing curved track extrapolation to get hits in road");
    }
469

470
    std::set<Identifier> addedIds;
471

472
473
474
475
476
    Muon::MdtPrepDataCollection::const_iterator pit = mdtPrdCol->begin();
    Muon::MdtPrepDataCollection::const_iterator pit_end = mdtPrdCol->end();
    for (; pit != pit_end; ++pit) {
        const Muon::MdtPrepData& mdtPrd = **pit;  // hit
        const Identifier& id = mdtPrd.identify();
477

478
479
        bool isFirst = isFirstProjection(id);
        Trk::MuonTrackSummary::ChamberHitSummary::Projection& proj = isFirst ? chamberHitSummary.m_first : chamberHitSummary.m_second;
480

481
        const Trk::Surface& surf = mdtPrd.detectorElement()->surface(id);
482

483
484
485
486
487
488
489
        const Trk::TrackParameters* exPars = 0;
        if (pars.associatedSurface() == surf) {
            exPars = &pars;
        } else {
            exPars = extrapolator->extrapolateDirectly(pars, surf, Trk::anyDirection, false, Trk::muon);
            if (!exPars) {
                if (isStraightLine) {
490
                    ATH_MSG_DEBUG(" Straight line propagation to prd " << m_idHelperSvc->toString(id) << " failed");
491
                } else {
492
                    ATH_MSG_DEBUG(" Curved track propagation to prd " << m_idHelperSvc->toString(id) << " failed");
493
494
495
                }
                continue;
            }
496
        }
497

498
499
        // use exPars to get distance to wire
        double distance = exPars->parameters()[Trk::locR];
500

501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
        // sometimes there is more than one hit in a tube,
        // which means there are two hits where the distance is the same but the tdc is different
        if (addedIds.count(id)) {
            ATH_MSG_DEBUG(" same tube hit, not adding to close hits in road");
        } else {
            // add all hits within the road width (defined in job options)
            if (std::abs(distance) < m_roadWidth) {
                ATH_MSG_VERBOSE("Hit ID  within road: " << m_idHelperSvc->toString(id) << " distance " << distance << " < " << m_roadWidth);
                ++proj.ncloseHits;
                addedIds.insert(id);
            } else {
                ATH_MSG_VERBOSE("Hit ID outside road: " << m_idHelperSvc->toString(id) << " distance " << distance
                                                        << " >= " << m_roadWidth);
            }
        }
        // to avoid double deleting when track is deleted, only delete
        // exPars when it's not the TrackParameters which was passed (pars)
        if (exPars != &pars) delete exPars;
519
520
    }

521
522
523
    // subtract the hits on the track in both projections:
    chamberHitSummary.m_first.ncloseHits -= chamberHitSummary.m_first.nhits;
    chamberHitSummary.m_second.ncloseHits -= chamberHitSummary.m_second.nhits;
524

525
526
527
528
529
    if (chamberHitSummary.m_first.ncloseHits < 0) {
        ATH_MSG_DEBUG("Number of hits in road < 0 in first projection: " << chamberHitSummary.m_first.ncloseHits
                                                                         << ", setting = 0. (nhits in first projection = "
                                                                         << chamberHitSummary.m_first.ncloseHits << ")");
        chamberHitSummary.m_first.ncloseHits = 0;
530
531
    }

532
533
534
535
536
537
538
    if (chamberHitSummary.m_second.ncloseHits < 0) {
        ATH_MSG_DEBUG("Number of hits in road < 0 in second projection: " << chamberHitSummary.m_second.ncloseHits
                                                                          << ", setting = 0. (nhits in second projection = "
                                                                          << chamberHitSummary.m_second.ncloseHits << ")");
        chamberHitSummary.m_second.ncloseHits = 0;
    }
}
539

540
541
542
543
bool Muon::MuonTrackSummaryHelperTool::isFirstProjection(const Identifier& id) const {
    if (!m_idHelperSvc->isMdt(id)) { return !m_idHelperSvc->measuresPhi(id); }
    return m_idHelperSvc->mdtIdHelper().multilayer(id) == 1;
}
544

545
546
const Muon::MdtPrepDataCollection* Muon::MuonTrackSummaryHelperTool::findMdtPrdCollection(const Identifier& chId) const {
    SG::ReadHandle<Muon::MdtPrepDataContainer> mdtPrdContainer(m_mdtKey);
547

548
549
550
    if (!mdtPrdContainer.isValid()) {
        ATH_MSG_WARNING("Cannot retrieve mdtPrepDataContainer " << m_mdtKey);
        return nullptr;
551
552
    }

553
554
555
    if (!mdtPrdContainer.isPresent()) {
        ATH_MSG_DEBUG("No MDT PRD container available");
        return nullptr;
556
557
    }

558
559
    IdentifierHash hash_id;
    m_idHelperSvc->mdtIdHelper().get_module_hash(chId, hash_id);
560

561
562
563
564
565
566
    auto coll = mdtPrdContainer->indexFindPtr(hash_id);
    if (coll == nullptr) {
        ATH_MSG_DEBUG(" MdtPrepDataCollection for:   " << m_idHelperSvc->toStringChamber(chId) << "  not found in container ");
        return nullptr;
    }
    return coll;
567
}