802 state.nscatterers = 0;
803 state.numberOfFlippedMdts = 0;
804 state.numberOfCleanedCompROTs = 0;
807 state.pullSum = ChamberPullInfo();
808 state.pullSumPhi = ChamberPullInfo();
809 state.pullSumTrigEta = ChamberPullInfo();
810 state.chambersToBeRemoved.clear();
811 state.chambersToBeRemovedPhi.clear();
812 state.pullSumPerChamber.clear();
813 state.pullSumPerChamberPhi.clear();
814 state.pullSumPerChamberEta.clear();
815 state.hitsPerChamber.clear();
816 state.outBoundsPerChamber.clear();
817 state.measInfo.clear();
818 state.largePullMeasurements.clear();
819 state.chamberLayerStatistics.clear();
820 state.stations.clear();
821 state.phiLayers.clear();
822 state.hasOfBoundsOutliers =
false;
823 state.hasVertexConstraint =
false;
824 state.hasSmall =
false;
825 state.hasLarge =
false;
827 state.nPseudoMeasurements = 0;
828 state.nPhiConstraints = 0;
833 if (!readHandle.isValid()) {
837 readHandle->getInitializedCache(fieldCache);
851 state.measInfo.reserve(
states->size());
853 std::set<int> rpcLayers;
854 std::set<int> tgcLayers;
856 double largestmdtpull = -999;
861 ATH_MSG_DEBUG(tsit->dumpType() <<
", parameters " << *tsit->trackParameters());
864 state.measInfo.emplace_back(tsit);
883 state.measInfo.emplace_back(tsit);
890 state.measInfo.emplace_back(tsit);
895 bool pseudo = !
id.is_valid();
896 if (pseudo) ++state.nPseudoMeasurements;
899 ATH_MSG_VERBOSE(
" TSOS is not a muon hit, position: r " <<
pars->position().perp() <<
" z " <<
pars->position().z());
903 state.measInfo.emplace_back(tsit);
909 state.hasVertexConstraint =
true;
913 bool measuresPhi = pseudo ? true :
m_idHelperSvc->measuresPhi(
id);
921 bool inBounds =
true;
923 double tol2 = 2 * tol1;
929 const MMClusterOnTrack* mmClusterOnTrack =
dynamic_cast<const MMClusterOnTrack*
>(meas);
934 }
else if (mmClusterOnTrack) {
936 inBounds = mmClusterOnTrack->detectorElement()->insideActiveBounds(
id, locPos, tol1, tol2);
947 ATH_MSG_DEBUG(
" calculation of residual/pull failed !!!!! ");
950 int pullSize = resPull->pull().size();
951 double residual = resPull->residual().front();
952 double pull = std::abs(resPull->pull().front());
955 if (!pseudo && pullSize != 1) {
964 double rTrackAbs = std::abs(rTrack);
965 double flippedResidual = isMDT ? rDrift + rTrack : 1e10;
966 double flippedPull = isMDT ? std::abs(flippedResidual /
error) : 1e10;
968 bool isNoise =
false;
971 bool flippedIsRecoverable =
973 double innerRadius = 14.6;
975 const MdtDriftCircleOnTrack* mdtdc =
dynamic_cast<const MdtDriftCircleOnTrack*
>(meas);
976 if (mdtdc) innerRadius = mdtdc->detectorElement()->innerTubeRadius();
979 bool isDelta = isOutlier && isMDT && rTrackAbs < innerRadius && rTrackAbs > std::abs(rDrift);
984 if (rTrackAbs > innerRadius) inBounds =
false;
985 }
else if (
pull > 10.) {
990 std::unique_ptr<MdtDriftCircleOnTrack> mdtRotFlipped;
991 std::unique_ptr<const CompetingMuonClustersOnTrack> updatedCompRot;
992 bool flipSign =
false;
994 const MdtDriftCircleOnTrack* mdtRot = isMDT ?
dynamic_cast<const MdtDriftCircleOnTrack*
>(meas) :
nullptr;
995 if (mdtRot && mdtRot->prepRawData() && mdtRot->prepRawData()->adc() <
m_adcCut) {
998 isRecoverable =
false;
999 flippedIsRecoverable =
false;
1002 if (flippedIsRecoverable) {
1007 if (!isNoise && flipSign) {
1009 mdtRotFlipped = std::make_unique<MdtDriftCircleOnTrack>(*mdtRot);
1012 double rDriftFlip = mdtRotFlipped->localParameters()[
Trk::locR];
1013 int signRot = rDrift < 0 ? -1 : 1;
1014 int signRotFlip = rDriftFlip < 0 ? -1 : 1;
1015 if (rDrift != 0. && signRot == signRotFlip) {
1016 ATH_MSG_WARNING(
" failed to flip sign of MDT " << rDrift <<
" flip " << rDriftFlip);
1019 ATH_MSG_WARNING(
" failed to dynamic_cast measurement with MDT Identifier to a MdtDriftCircleOnTrack");
1024 largestmdtpull =
pull;
1038 rpcLayers.insert(
layer);
1047 if (stName[1] ==
'1')
1049 else if (stName[1] ==
'2')
1051 else if (stName[1] ==
'3')
1063 tgcLayers.insert(
layer);
1069 const CompetingMuonClustersOnTrack* crot = (measuresPhi && !isMDT &&
m_idHelperSvc->isRpc(
id))
1070 ?
dynamic_cast<const CompetingMuonClustersOnTrack*
>(meas)
1073 ATH_MSG_DEBUG(
" CompetingMuonClustersOnTrack with rots " << crot->numberOfContainedROTs());
1077 double absminres = 0.;
1078 double absmaxres = 0.;
1079 for (
unsigned int i = 0;
i < crot->numberOfContainedROTs(); ++
i) {
1080 const MuonClusterOnTrack* cluster = &crot->rioOnTrack(
i);
1081 if (!cluster)
continue;
1087 minpos = cluster->localParameters()[
Trk::locX];
1090 }
else if (absres < absminres) {
1093 minpos = cluster->localParameters()[
Trk::locX];
1094 }
else if (absres > absmaxres) {
1102 ATH_MSG_DEBUG(
" residuals: min " << minres <<
" max " << maxres <<
" diff " << maxres - minres);
1103 bool splitCompRot =
false;
1104 if (std::abs(maxres - minres) > 100 && absmaxres - absminres > 20 && crot->numberOfContainedROTs() < 20) {
1106 splitCompRot =
true;
1109 std::list<const Trk::PrepRawData*> prdList;
1111 for (
unsigned int i = 0;
i < crot->numberOfContainedROTs(); ++
i) {
1112 const MuonClusterOnTrack* cluster = &crot->rioOnTrack(
i);
1113 if (!cluster)
continue;
1118 << cluster->localParameters()[
Trk::locX] <<
" pars "
1120 prdList.push_back(cluster->prepRawData());
1123 if (prdList.empty()) {
1124 ATH_MSG_WARNING(
"No clusters selected during comprot cleaning, keeping old cluster");
1127 ++state.numberOfCleanedCompROTs;
1134 state.measInfo.emplace_back(
id, chId, chIndex, inBounds,
residual,
pull, tsit, meas,
pars, std::move(resPull),
nullptr);
1135 MCTBCleaningInfo&
info = state.measInfo.back();
1136 if (flipSign) {
info.flippedMdt = std::move(mdtRotFlipped); }
1137 info.isNoise = isNoise;
1138 if (updatedCompRot) {
1140 info.cleanedCompROT = std::move(updatedCompRot);
1142 std::unique_ptr<Trk::TrackParameters> exPars =
m_extrapolator->extrapolate(ctx,
1143 *
pars,
info.cleanedCompROT->associatedSurface(),
1146 ATH_MSG_WARNING(
"Update of comp rot parameters failed, keeping old ones");
1147 info.cleanedCompROT.reset();
1149 info.pars = exPars.get();
1150 state.parsToBeDeleted.emplace_back(std::move(exPars));
1155 ChamberLayerStatistics& chamberStatistics = state.chamberLayerStatistics[chIndex];
1156 chamberStatistics.chIndex = chIndex;
1161 ATH_MSG_DEBUG(
" r_drift " << rDrift <<
" r_trk " << rTrack <<
" dr " << std::setw(7)
1168 if (flippedIsRecoverable)
1170 else if (isRecoverable)
1177 state.hasOfBoundsOutliers =
true;
1178 ++chamberStatistics.noutBounds;
1179 }
else if (flippedIsRecoverable || isRecoverable)
1180 ++chamberStatistics.nrecoverableOutliers;
1182 ++chamberStatistics.ndeltas;
1184 ++chamberStatistics.noutliers;
1186 if (!pseudo) ++state.noutliers;
1192 ++chamberStatistics.nhits;
1193 if (!pseudo) ++state.nhits;
1195 if (flipSign) ++state.numberOfFlippedMdts;
1197 std::string idString = pseudo ?
" Pseudomeasurement " :
m_idHelperSvc->toString(
id);
1198 std::string boundCheck = inBounds ?
"inBounds" :
"outBounds";
1205 ATH_MSG_VERBOSE(
" r_drift " << std::setw(8) << rDrift <<
" r_trk " << std::setw(7) << rTrack <<
" dr " << std::setw(7)
1212 const RpcClusterOnTrack* rpc =
dynamic_cast<const RpcClusterOnTrack*
>(meas);
1214 const CompetingMuonClustersOnTrack* crot =
dynamic_cast<const CompetingMuonClustersOnTrack*
>(meas);
1215 if (crot) { rpc =
dynamic_cast<const RpcClusterOnTrack*
>(crot->containedROTs().front()); }
1217 if (rpc)
ATH_MSG_DEBUG(
" time " << rpc->prepRawData()->time() - rpc->globalPosition().mag() / 300.);
1222 if (!isNoise && isOutlier) { state.largePullPseudoMeasurements.insert(meas); }
1227 EtaPhiHits&
hits = state.hitsPerChamber[
info.chId];
1235 EtaPhiHits& outHits = state.outBoundsPerChamber[
info.chId];
1247 ChamberPullInfo& pullInfoHits = measuresPhi ? state.pullSumPhi : state.pullSum;
1248 pullInfoHits.pullSum +=
pull;
1249 ++pullInfoHits.nhits;
1250 if (
pull > pullInfoHits.maxPull) pullInfoHits.maxPull =
pull;
1252 ChamberPullInfo& pullInfoCh = measuresPhi ? state.pullSumPerChamberPhi[
info.chId] : state.pullSumPerChamber[
info.chId];
1253 pullInfoCh.pullSum +=
pull;
1255 if (
pull > pullInfoCh.maxPull) pullInfoCh.maxPull =
pull;
1257 ChamberPullInfo& pullInfoTrig = measuresPhi ? state.pullSumPhi : state.pullSumTrigEta;
1258 if (!measuresPhi) state.pullSumPerChamberEta[
info.chId];
1259 pullInfoTrig.pullSum +=
pull;
1260 ++pullInfoTrig.nhits;
1261 if (
pull > pullInfoTrig.maxPull) pullInfoTrig.maxPull =
pull;
1267 if (mdtmeas) state.largePullMeasurements.insert(mdtmeas);
1270 unsigned int nphiLayers = rpcLayers.size() + tgcLayers.size();
1271 if (state.hasVertexConstraint) { nphiLayers += 1; }
1272 if (state.nIdHits > 0) { nphiLayers += 2; }
1274 if ((nphiLayers > 2) && state.numberOfCleanedCompROTs > 0) {
1275 ATH_MSG_DEBUG(
" Sufficient constraints to clean competing ROTs: trigger layers " << nphiLayers);
1278 state.numberOfCleanedCompROTs = 0;
1282 if (!
m_use_slFit && state.slFit && (state.hasVertexConstraint || state.nIdHits > 0) && fieldCache.
solenoidOn()) {
1283 state.slFit =
false;
1286 if (state.hasVertexConstraint || state.nIdHits) {
1287 if (state.hasVertexConstraint)
ATH_MSG_DEBUG(
" Track has vertex contraint:");
1288 if (state.nIdHits > 0)
ATH_MSG_DEBUG(
" Track has ID Hits: " << state.nIdHits);
1303 std::set<Identifier> chambersToBeRemoved, chambersToBeRemovedPhi, goodEtaLayers, goodPhiLayers;
1315 for (; chit != chit_end; ++chit) {
1316 if (!chit->first.is_valid())
continue;
1322 state.stations.insert(stIndex);
1325 EtaPhiHits& nhits = chit->second;
1326 EtaPhiHits& noutBounds = state.outBoundsPerChamber[chit->first];
1328 << noutBounds.neta <<
" phi hits " << nhits.nphi <<
" outBounds " << noutBounds.nphi);
1329 if (nhits.neta != 0 && nhits.neta == noutBounds.neta) {
1332 double fakePull = 1000. * nhits.neta;
1333 if (!chambersToBeRemoved.count(chit->first)) {
1335 state.chambersToBeRemoved.emplace_back(fakePull, chit->first);
1336 chambersToBeRemoved.insert(chit->first);
1339 PullChIt pIt = state.chambersToBeRemoved.begin();
1340 PullChIt pIt_end = state.chambersToBeRemoved.end();
1341 for (; pIt != pIt_end; ++pIt) {
1342 if (pIt->second == chit->first && pIt->first < fakePull) {
1343 pIt->first = fakePull;
1349 if (noutBounds.neta != 0)
ATH_MSG_DEBUG(
" --> Some eta hits out of bounds ");
1351 if (isPrec && nhits.neta > 0) {
1353 state.hasSmall =
true;
1355 state.hasLarge =
true;
1358 if (nhits.nphi != 0 && nhits.nphi == noutBounds.nphi) {
1361 double fakePull = 1000. * nhits.nphi;
1362 if (!chambersToBeRemovedPhi.count(chit->first)) {
1364 state.chambersToBeRemovedPhi.emplace_back(fakePull, chit->first);
1365 chambersToBeRemovedPhi.insert(chit->first);
1368 PullChIt pIt = state.chambersToBeRemovedPhi.begin();
1369 PullChIt pIt_end = state.chambersToBeRemovedPhi.end();
1370 for (; pIt != pIt_end; ++pIt) {
1371 if (pIt->second == chit->first) {
1372 pIt->first =
std::max(pIt->first, fakePull);
1378 if (noutBounds.nphi != 0)
ATH_MSG_DEBUG(
" --> Some phi hits out of bounds ");
1379 if (nhits.nphi > 0) {
1386 state.nPhiConstraints = state.phiLayers.size() + state.nPseudoMeasurements;
1387 if (state.hasLarge && state.hasSmall) ++state.nPhiConstraints;
1388 if (state.stations.size() == 1 && nphiLayers > 1) ++state.nPhiConstraints;
1390 double pull_precisionhits = -1.;
1391 if (state.pullSum.nhits != 0) { pull_precisionhits = state.pullSum.pullSum / state.pullSum.nhits; }
1392 double pull_triggerhits = -1.;
1393 if (state.pullSumTrigEta.nhits != 0) { pull_triggerhits = state.pullSumTrigEta.pullSum / state.pullSumTrigEta.nhits; }
1394 double pull_phihits = -1.;
1395 if (state.pullSumPhi.nhits != 0) { pull_phihits = state.pullSumPhi.pullSum / state.pullSumPhi.nhits; }
1397 << pull_precisionhits <<
" trigger eta " << pull_triggerhits <<
" phi " << pull_phihits <<
" measurements "
1398 << state.measInfo.size() << std::endl
1399 <<
" precision chambers " << state.pullSumPerChamber.size() <<
" hits with large pull "
1400 << state.largePullMeasurements.size() <<
" pseudos with large pull " << state.largePullPseudoMeasurements.size()
1401 <<
" chambers to be removed " << state.chambersToBeRemoved.size() <<
" phi " << state.chambersToBeRemovedPhi.size()
1402 <<
" phi lay " << state.phiLayers.size() <<
" phi constr " << state.nPhiConstraints);
1404 if (state.pullSumPhi.nhits != 0) {
1405 double phiPull = state.pullSumPhi.pullSum / state.pullSumPhi.nhits;
1406 if (state.stations.size() == 1 && phiPull > 5.) {
ATH_MSG_DEBUG(
" Single station track with large phi pulls!! "); }