1008 {
1010
1011
1012
1014 std::vector<Tracklet> myTracks = tracks;
1015 tracks.clear();
1016 for (const Tracklet &track : myTracks) {
1017 Identifier id1 =
track.getML1seg().mdtHitsOnTrack().at(0)->identify();
1018 Identifier
id2 =
track.getML2seg().mdtHitsOnTrack().at(0)->identify();
1019 int nLayerML1 =
m_idHelperSvc->mdtIdHelper().tubeLayerMax(id1);
1021 double ratio = (
double)(
track.mdtHitsOnTrack().size()) / (nLayerML1 + nLayerML2);
1022 if (ratio > 0.75) tracks.push_back(track);
1023 }
1024 }
1025
1026 std::vector<Tracklet> UniqueTracks;
1027 std::vector<unsigned int> AmbigTrks;
1028 for (unsigned int tk1 = 0; tk1 < tracks.size(); ++tk1) {
1029 int nShared = 0;
1030
1031 bool isResolved = false;
1032 for (unsigned int AmbigTrksIdx : AmbigTrks) {
1033 if (tk1 == AmbigTrksIdx) {
1034 isResolved = true;
1035 break;
1036 }
1037 }
1038 if (isResolved) continue;
1039 std::vector<Tracklet> AmbigTracks;
1040 AmbigTracks.push_back(tracks.at(tk1));
1041
1042 double Trk1ML1R = tracks.at(tk1).getML1seg().globalPosition().perp();
1043 double Trk1ML1Z = tracks.at(tk1).getML1seg().globalPosition().z();
1044 double Trk1ML2R = tracks.at(tk1).getML2seg().globalPosition().perp();
1045 double Trk1ML2Z = tracks.at(tk1).getML2seg().globalPosition().z();
1046
1047 Identifier tk1ID = tracks.at(tk1).muonIdentifier();
1048 bool tk1_isBarrel =
m_idHelperSvc->mdtIdHelper().isBarrel(tk1ID);
1049 bool tk1_isEndcap =
m_idHelperSvc->mdtIdHelper().isEndcap(tk1ID);
1050
1051
1052 for (unsigned int tk2 = (tk1 + 1); tk2 < tracks.size(); ++tk2) {
1053 if (tracks.at(tk1).mdtChamber() == tracks.at(tk2).mdtChamber() && tracks.at(tk1).mdtChPhi() == tracks.at(tk2).mdtChPhi() &&
1054 (tracks.at(tk1).mdtChEta()) * (tracks.at(tk2).mdtChEta()) > 0) {
1055
1056 for (unsigned int AmbigTrksIdx : AmbigTrks) {
1057 if (tk2 == AmbigTrksIdx) {
1058 isResolved = true;
1059 break;
1060 }
1061 }
1062 if (isResolved) continue;
1063
1064 double Trk2ML1R = tracks.at(tk2).getML1seg().globalPosition().perp();
1065 double Trk2ML1Z = tracks.at(tk2).getML1seg().globalPosition().z();
1066 double Trk2ML2R = tracks.at(tk2).getML2seg().globalPosition().perp();
1067 double Trk2ML2Z = tracks.at(tk2).getML2seg().globalPosition().z();
1068
1069
1070 double DistML1(1000), DistML2(1000);
1071 if (tk1_isBarrel) {
1072 DistML1 = std::abs(Trk1ML1Z - Trk2ML1Z);
1073 DistML2 = std::abs(Trk1ML2Z - Trk2ML2Z);
1074 } else if (tk1_isEndcap) {
1075 DistML1 = std::abs(Trk1ML1R - Trk2ML1R);
1076 DistML2 = std::abs(Trk1ML2R - Trk2ML2R);
1077 }
1078 if (DistML1 < 40 || DistML2 < 40) {
1079
1080 std::vector<const Muon::MdtPrepData*> mdts1 = tracks.at(tk1).mdtHitsOnTrack();
1081 std::vector<const Muon::MdtPrepData*> mdts2 = tracks.at(tk2).mdtHitsOnTrack();
1082 nShared = 0;
1083 for (const Muon::MdtPrepData *mdt1 : mdts1) {
1084 for (const Muon::MdtPrepData *mdt2 : mdts2) {
1085 if (mdt1->identify() == mdt2->identify()) {
1086 ++nShared;
1087 break;
1088 }
1089 }
1090 }
1091
1092 if (nShared <= 1) continue;
1093
1094 AmbigTracks.push_back(tracks.at(tk2));
1095 AmbigTrks.push_back(tk2);
1096 }
1097 }
1098 }
1099
1100 if (AmbigTracks.size() == 1) {
1101 UniqueTracks.push_back(tracks.at(tk1));
1102 continue;
1103 }
1104
1105
1106 if (tk1_isBarrel) {
1107 bool hasMomentum = tracks.at(tk1).charge() != 0;
1108 double aveX(0), aveY(0), aveZ(0), aveAlpha(0);
1109 double aveP(0), nAmbigP(0), TrkCharge(tracks.at(tk1).charge());
1110 bool allSameSign(true);
1111
1112 for (const Tracklet &AmbigTrack : AmbigTracks) {
1113 if (!hasMomentum) {
1114 aveX += AmbigTrack.globalPosition().x();
1115 aveY += AmbigTrack.globalPosition().y();
1116 aveZ += AmbigTrack.globalPosition().z();
1117 aveAlpha += AmbigTrack.getML1seg().alpha();
1118 } else {
1119
1120 if (std::abs(AmbigTrack.charge() - TrkCharge) > 0.1) allSameSign = false;
1121
1122 aveP += AmbigTrack.momentum().mag();
1123 ++nAmbigP;
1124 aveAlpha += AmbigTrack.alpha();
1125 aveX += AmbigTrack.globalPosition().x();
1126 aveY += AmbigTrack.globalPosition().y();
1127 aveZ += AmbigTrack.globalPosition().z();
1128 }
1129 }
1130 if (!hasMomentum) {
1131 aveX = aveX / (
double)AmbigTracks.size();
1132 aveY = aveY / (
double)AmbigTracks.size();
1133 aveZ = aveZ / (
double)AmbigTracks.size();
1135 aveAlpha = aveAlpha / (
double)AmbigTracks.size();
1136 double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1137 double rErr = tracks.at(tk1).getML1seg().rError();
1138 double zErr = tracks.at(tk1).getML1seg().zError();
1139
1140 TrackletSegment aveSegML1(
m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1141 double pT =
m_maxpTot * std::sin(aveSegML1.alpha());
1142 double pz =
m_maxpTot * std::cos(aveSegML1.alpha());
1144 pT * std::sin(aveSegML1.globalPosition().phi()),
1145 pz);
1147 matrix.setIdentity();
1148 matrix(0, 0) = std::
pow(tracks.at(tk1).getML1seg().rError(),2);
1149 matrix(1, 1) = std::
pow(tracks.at(tk1).getML1seg().zError(),2);
1150 matrix(2, 2) = std::
pow(0.0000001,2);
1151 matrix(3, 3) = std::
pow(tracks.at(tk1).getML1seg().alphaError(),2);
1153 Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1154 UniqueTracks.push_back(aveTrack);
1155 } else
if (allSameSign) {
1156 aveP = aveP / nAmbigP;
1157 double pT = aveP * std::sin(tracks.at(tk1).getML1seg().alpha());
1158 double pz = aveP * std::cos(tracks.at(tk1).getML1seg().alpha());
1160 pT * std::sin(tracks.at(tk1).globalPosition().phi()),
1161 pz);
1162 Tracklet MyTrack = tracks.at(tk1);
1164 MyTrack.
charge(tracks.at(tk1).charge());
1165 UniqueTracks.push_back(MyTrack);
1166 } else {
1167 aveX = aveX / (
double)AmbigTracks.size();
1168 aveY = aveY / (
double)AmbigTracks.size();
1169 aveZ = aveZ / (
double)AmbigTracks.size();
1171 aveAlpha = aveAlpha / (
double)AmbigTracks.size();
1172 double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1173 double rErr = tracks.at(tk1).getML1seg().rError();
1174 double zErr = tracks.at(tk1).getML1seg().zError();
1175
1176 TrackletSegment aveSegML1(
m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1177 double pT =
m_maxpTot * std::sin(aveSegML1.alpha());
1178 double pz =
m_maxpTot * std::cos(aveSegML1.alpha());
1180 pT * std::sin(aveSegML1.globalPosition().phi()),
1181 pz);
1183 matrix.setIdentity();
1184 matrix(0, 0) = std::
pow(tracks.at(tk1).getML1seg().rError(),2);
1185 matrix(1, 1) = std::
pow(tracks.at(tk1).getML1seg().zError(),2);
1186 matrix(2, 2) = std::
pow(0.0000001,2);
1187 matrix(3, 3) = std::
pow(tracks.at(tk1).getML1seg().alphaError(),2);
1189 Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1190 UniqueTracks.push_back(aveTrack);
1191 }
1192 }
1193
1194
1195 else
if (tk1_isEndcap) {
1196 std::vector<const Muon::MdtPrepData*> AllMdts;
1197 for (Tracklet const &AmbigTrack : AmbigTracks) {
1198 std::vector<const Muon::MdtPrepData*> mdts = AmbigTrack.mdtHitsOnTrack();
1199 std::vector<const Muon::MdtPrepData*> tmpAllMdt = AllMdts;
1200 for (const Muon::MdtPrepData *mdt : mdts) {
1201 bool isNewHit = true;
1202 for (const Muon::MdtPrepData *tmpmdt : tmpAllMdt) {
1203 if (mdt->identify() == tmpmdt->identify()) {
1204 isNewHit = false;
1205 break;
1206 }
1207 }
1208 if (isNewHit) AllMdts.push_back(mdt);
1209 }
1210 }
1211
1213 if (!MyECsegs.empty()) {
1214 TrackletSegment ECseg = MyECsegs.at(0);
1220 pz);
1222 matrix.setIdentity();
1223 matrix(0, 0) = std::
pow(ECseg.rError(),2);
1224 matrix(1, 1) = std::
pow(ECseg.zError(),2);
1225 matrix(2, 2) = std::
pow(0.0000001,2);
1226 matrix(3, 3) = std::
pow(ECseg.alphaError(),2);
1228 Tracklet MyCombTrack(MyECsegs.at(0), ECseg, momentum, matrix, 0);
1229 UniqueTracks.push_back(MyCombTrack);
1230 } else
1231 UniqueTracks.push_back(tracks.at(tk1));
1232 }
1233
1234 }
1235
1236 return UniqueTracks;
1237 }
void momentum(const Amg::Vector3D &p)
void charge(double charge)