25 using MuonSegmentVec = std::vector<std::unique_ptr<Muon::MuonSegment>>;
47 std::stringstream sstr{};
52 constexpr
double minEtaNSW = 1.2;
53 constexpr
double maxEtaNSW = 2.8;
105 m_dir{
cl->detectorElement()->transform(
cl->identify()).linear() * Amg::Vector3D::UnitY()} {}
108 const std::array<double,2>& lengths) :
110 m_pos{seed[0].pos() + lengths[0] * seed[0].dir()} {
113 m_dir = un_dir.unit()*(un_dir.z() * seed[0].pos().z() > 0 ? 1 : -1);
119 for (
const SeedMeasurement&
cl : seed) {
135 m_parent{
parent},
m_pos{seg.globalPosition()},
m_dir{seg.globalDirection()} {
152 const double BdotD = B.dot(D);
155 if (std::abs(
divisor) < std::numeric_limits<double>::epsilon())
return diff.mag();
158 return (
beta * B - delta * D -
diff).mag();
161 if (!
size() || !meas)
return false;
162 if (
find(meas))
return true;
178 if (!prd)
return false;
180 if (!design)
return false;
182 const double dist = std::hypot(padDist.x(), padDist.y());
198 unsigned int only_th{0}, only_oth{0};
206 if (!
other.m_measurements[
m]) ++only_th;
207 if (!
other.m_phiMeasurements[
m]) ++only_th;
215 meas.reserve(
size());
222 return std::abs(a->globalPosition().z()) < std::abs(b->globalPosition().z());
229 return insert(std::move(meas));
235 if (!seed || meas.
distance() < seed.distance()) {
246 return m_calibClust.emplace(std::move(new_clust)).first->get();
258 declareInterface<IMuonNSWSegmentFinderTool>(
this);
273 ATH_MSG_FATAL(
"The muon should come from the IP && horizontally from the Calo..."<<
274 "We need to place a very good basball player at the calo exit to deflect the muon horizontally.");
275 return StatusCode::FAILURE;
278 return StatusCode::SUCCESS;
285 std::vector<const Muon::MuonClusterOnTrack*> muonClusters{};
286 muonClusters.reserve(cache.
inputClust.size());
288 [](
const std::unique_ptr<const MuonClusterOnTrack>&
cl){return cl.get();});
289 ATH_MSG_DEBUG(
"Entering MuonNSWSegmentFinderTool with " << muonClusters.size() <<
" clusters to be fit");
292 MuonSegmentVec out_segments{};
295 ATH_MSG_VERBOSE(
"Found " << stereoSegs.size() <<
" MMG stereo seeded segments");
296 out_segments.insert(out_segments.end(), std::make_move_iterator(stereoSegs.begin()),
297 std::make_move_iterator(stereoSegs.end()));
301 auto dump_output = [&]() {
303 for (std::unique_ptr<Muon::MuonSegment>& seg : out_segments) {
312 std::vector<const Muon::MuonClusterOnTrack*> clustPostStereo{};
314 std::array<std::set<Identifier>, 16> masked_segs{};
315 for (std::unique_ptr<Muon::MuonSegment>& seg : out_segments) {
323 if (!out_segments.empty()) {
324 clustPostStereo.reserve(muonClusters.size());
332 const std::vector<const Muon::MuonClusterOnTrack*>& segmentInput = !out_segments.empty() ? clustPostStereo : muonClusters;
337 ATH_MSG_VERBOSE(
"Found " << etaSegs.size() <<
" stgc seeded eta segments");
338 MuonSegmentVec precSegs =
find3DSegments(ctx, segmentInput, etaSegs);
340 out_segments.insert(out_segments.end(),
341 std::make_move_iterator(precSegs.begin()),
342 std::make_move_iterator(precSegs.end()));
352 for (std::unique_ptr<Muon::MuonSegment>& seg : out_segments) {
353 std::vector<const MuonClusterOnTrack*> seg_hits{};
354 seg_hits.reserve(seg->containedMeasurements().size());
357 if (clus) seg_hits.push_back(clus);
360 for (
int iWedge{1}; iWedge<=4 ; ++iWedge) {
363 if ( quad_hits.size () < 2 || ((iWedge == 2 || iWedge == 3) && quad_hits.size() < 4))
continue;
364 std::vector<const Trk::MeasurementBase*> fit_meas{};
365 std::unique_ptr<Trk::PseudoMeasurementOnTrack> pseudoVtx{
ipConstraint(ctx)};
366 if (pseudoVtx) fit_meas.push_back(pseudoVtx.get());
367 std::copy(quad_hits.begin(), quad_hits.end(), std::back_inserter(fit_meas));
369 const Trk::Surface& surf = quad_hits.front()->associatedSurface();
374 if (perpos.dot(gdir_seg) < 0) gdir_seg *= -1;
377 std::unique_ptr<Trk::Track> segtrack =
fit(ctx, fit_meas, startpar);
378 if (!segtrack)
continue;
383 <<
"position: "<<
to_string(seg->globalPosition())<< std::endl
384 <<
"direction: "<<
to_string(seg->globalDirection())<< std::endl
385 <<
m_printer->print(seg->containedMeasurements()));
387 cache.
quadSegs.emplace_back(std::move(seg));
394 const std::vector<const Muon::MuonClusterOnTrack*>& allClusts,
395 int singleWedge)
const {
407 for (
MeasVec& hitsInLayer : orderedClust) hitsInLayer =
vetoBursts(std::move(hitsInLayer));
408 orderedClust.erase(std::remove_if(orderedClust.begin(), orderedClust.end(),
409 [](
const MeasVec&
vec) { return vec.empty(); }),
411 if (orderedClust.empty())
return {};
415 ATH_MSG_DEBUG(
"Retrieved " << seeds.size() <<
" seeds in the MMStereoAlg after the ambiguity resolution" );
417 if (seeds.empty())
return {};
423 std::vector<const Trk::MeasurementBase*> fit_meas{};
425 std::unique_ptr<Trk::PseudoMeasurementOnTrack> pseudoVtx{
ipConstraint(ctx)};
426 if (pseudoVtx) fit_meas.push_back(pseudoVtx.get());
428 std::copy(calib_clust.begin(), calib_clust.end(), std::back_inserter(fit_meas));
430 const Trk::Surface& surf = calib_clust.front()->associatedSurface();
436 if (perpos.dot(gdir_seg) < 0) gdir_seg *= -1;
441 std::unique_ptr<Trk::Track> segtrack =
fit(ctx, fit_meas, startpar);
442 if (segtrack) trackSegs.push_back(std::move(segtrack));
448 const std::vector<const Trk::MeasurementBase*>& fit_meas,
451 to_string(perigee.
momentum())<<
". Contained measurements in candidate: " << std::endl
459 std::unique_ptr<Trk::Track> cleanedTrack =
m_trackCleaner->clean(*segtrack, ctx);
479 const std::vector<const Muon::MuonClusterOnTrack*>& muonClusters,
480 int singleWedge)
const {
493 ATH_MSG_DEBUG(
"Cleaning eta clusters in stgc seeded 2D segments ");
500 if (orderedClusters.size() < 4)
return {};
510 if (seed.size() < 4){
511 ATH_MSG_VERBOSE(__func__<<
" :"<<__LINE__<<
" - Seed size: " << seed.size() <<
" is below the cut (4)");
515 etaHitVec = seed.measurements();
533 if (perpos.dot(gdir_seg) < 0) gdir_seg *= -1;
534 const auto startpar =
Trk::Perigee(perpos, gdir_seg, 0, perpos);
535 ATH_MSG_VERBOSE(
" start parameter " << perpos <<
" pp " << startpar.position() <<
" gd " << gdir_seg.unit() <<
" pp "
536 << startpar.momentum().unit());
539 hitsToTrack(ctx, etaHitVec, phiHitVec, startpar, segTrkColl);
556 MuonSegmentVec segments{};
557 std::unique_ptr<const TrackCollection> resolvedTracks(
m_ambiTool->process(&segTrkColl));
558 ATH_MSG_DEBUG(
"Resolved track candidates: old size " << segTrkColl.
size() <<
" new size " << resolvedTracks->
size());
561 for (
const Trk::Track* trk : *resolvedTracks) {
562 const auto* measurements = trk->measurementsOnTrack();
563 const bool has_eta = std::find_if(measurements->begin(), measurements->end(),
565 Identifier id = m_edmHelperSvc->getIdentifier(*meas);
566 return id.is_valid() && !m_idHelperSvc->measuresPhi(id);
567 }) != trk->measurementsOnTrack()->end();
568 if (!has_eta)
continue;
573 segments.emplace_back(std::move(seg));
582 const std::vector<const Muon::MuonClusterOnTrack*>& muonClusters,
583 MuonSegmentVec& etaSegs,
584 int singleWedge)
const {
585 MuonSegmentVec segments{};
588 ATH_MSG_DEBUG(
"Cleaning phi clusters in stgc seeded 3D segments ");
590 ATH_MSG_DEBUG(
"After hit cleaning, there are " << phiClusters.size() <<
" phi clusters to be fit");
596 if (orderedWireClusters.size() + orderedPadClusters.size() < 2) {
597 ATH_MSG_DEBUG(
"Not enough phi hits present, cannot perform the 3D fit!");
598 segments.insert(segments.end(), std::make_move_iterator(etaSegs.begin()), std::make_move_iterator(etaSegs.end()));
603 ATH_MSG_DEBUG(
"Cleaning eta clusters in stgc seeded 3D segments ");
608 bool triedWireSeed{
false};
609 std::vector<NSWSeed> seeds_WiresSTGC;
612 for (std::unique_ptr<Muon::MuonSegment>& etaSeg : etaSegs) {
617 std::vector<NSWSeed> seeds;
620 if (std::abs(etaSeg->globalPosition().eta()) < 2.4) {
621 if (!triedWireSeed) {
623 triedWireSeed =
true;
627 if (!seeds_WiresSTGC.empty()) {
628 seeds = seeds_WiresSTGC;
634 if (seeds.empty() ) {
639 std::back_inserter(orderedPadClustersForSeeding),
649 double etaSegLocX = etaSeg->localParameters()[
Trk::locX];
650 double etaSegLocXZ = etaSeg->localDirection().angleXZ();
670 if (perpos.dot(gdir_seg) < 0) gdir_seg *= -1;
671 const Trk::Perigee startpar{perpos, gdir_seg, 0, perpos};
673 NSWSeed seed3D{
this, perpos, gdir_seg};
678 if (nPhiHits < 2)
continue;
680 MeasVec phiHitVec = seed3D.measurements();
685 if (
hitsToTrack(ctx, etaHitsCalibrated, phiHitVec, startpar, segTrkColl)) {
687 ATH_MSG_VERBOSE(
"Segment successfully fitted for wedge "<<singleWedge<<std::endl<<
694 if (!is3Dseg) { segments.push_back(std::move(etaSeg)); }
697 segments.insert(segments.end(), std::make_move_iterator(new_segs.begin()), std::make_move_iterator(new_segs.end()));
705 covVtx.setIdentity();
706 covVtx = errVtx * errVtx * covVtx;
722 covVtx(0, 0) = errVtx * errVtx;
747 std::vector<const Trk::MeasurementBase*> vecFitPts;
748 unsigned int nHitsEta = etaHitVec.size();
749 unsigned int nHitsPhi = phiHitVec.size();
752 std::unique_ptr<Trk::PseudoMeasurementOnTrack> pseudoVtx{
ipConstraint(ctx)},
754 pseudoPhi1{
nullptr}, pseudoPhi2{
nullptr};
756 if (pseudoVtx) { vecFitPts.push_back(pseudoVtx.get()); }
757 if (pseudoVtxCalo) { vecFitPts.push_back(pseudoVtxCalo.get()); }
765 cov(0, 0) = pseudoPrecision;
767 pseudoPhi1 = std::make_unique<Trk::PseudoMeasurementOnTrack>(
Trk::LocalParameters(loc_pseudopars),
769 etaHitVec.front()->associatedSurface());
770 pseudoPhi2 = std::make_unique<Trk::PseudoMeasurementOnTrack>(
Trk::LocalParameters(loc_pseudopars),
772 etaHitVec.back()->associatedSurface());
775 vecFitPts.push_back(pseudoPhi1.get());
776 std::copy(etaHitVec.begin(), etaHitVec.end(), std::back_inserter(vecFitPts));
777 vecFitPts.push_back(pseudoPhi2.get());
778 ATH_MSG_DEBUG(
"Fitting a 2D-segment track with " << nHitsEta <<
" Eta hits");
782 std::merge(phiHitVec.begin(), phiHitVec.end(), etaHitVec.begin(), etaHitVec.end(), std::back_inserter(vecFitPts),
784 double z1 = std::abs(c1->detectorElement()->center(c1->identify()).z());
785 double z2 = std::abs(c2->detectorElement()->center(c2->identify()).z());
788 ATH_MSG_DEBUG(
"Fitting a 3D-segment track with " << nHitsEta <<
" Eta hits and " << nHitsPhi <<
" Phi hits");
792 std::unique_ptr<Trk::Track> segtrack =
fit(ctx, vecFitPts, startpar);
793 if (!segtrack)
return false;
794 segTrkColl.
push_back(std::move(segtrack));
800 const std::vector<const Muon::MuonClusterOnTrack*>& muonClusters,
int hit_sel,
int singleWedge )
const {
804 clusters.reserve(muonClusters.size());
806 if (!cluster)
continue;
807 if (singleWedge && singleWedge !=
wedgeNumber(cluster))
continue;
826 std::array<std::set<Identifier>,16> used_hits{};
836 const int channelType =
m_idHelperSvc->stgcIdHelper().channelType(
id);
841 std::set<Identifier>& lay_hits = used_hits[iorder];
842 if (lay_hits.count(
id))
continue;
844 orderedClusters[iorder].emplace_back(hit);
846 if (
nBad)
ATH_MSG_WARNING(
"Unable to classify " <<
nBad <<
" clusters by their layer since they are neither MM nor sTGC");
849 orderedClusters.erase(std::remove_if(orderedClusters.begin(), orderedClusters.end(),
850 [](
const MeasVec&
vec) { return vec.empty(); }),
851 orderedClusters.end());
854 for(
MeasVec& lays: orderedClusters){
856 return channel(a) < channel(b);
864 return orderedClusters;
870 std::vector<NSWSeed> seeds;
874 if (orderedClusters.size() < 4)
return seeds;
880 int seedingLayersL{0};
881 for (
unsigned int ilayerL{0}; (ilayerL < orderedClusters.size() && seedingLayersL <
m_nOfSeedLayers); ++ilayerL) {
882 bool usedLayerL{
false};
885 if (usePhi !=
m_idHelperSvc->measuresPhi(hitL->identify()))
continue;
890 int seedingLayersR{0};
891 for (
unsigned int ilayerR = orderedClusters.size() - 1; (ilayerR > ilayerL && seedingLayersR <
m_nOfSeedLayers);
893 bool usedLayerR{
false};
895 if (usePhi !=
m_idHelperSvc->measuresPhi(hitR->identify()))
continue;
899 const double eta = seed.dir().perp() > std::numeric_limits<float>::epsilon() ? std::abs(seed.dir().eta()): FLT_MAX;
900 if (eta < minEtaNSW || eta > maxEtaNSW) {
905 if(std::abs(std::abs(
m_idHelperSvc->stationEta(hitL->identify())) - std::abs(
m_idHelperSvc->stationEta(hitR->identify()))) > 1) {
909 Trk::Intersection intersect = hitR->detectorElement()->surface(hitR->identify()).straightLineIntersection(hitL->globalPosition(),hitL->globalPosition().unit(),
false,
false);
911 hitR->detectorElement()->surface(hitR->identify()).globalToLocal(
intersect.position,
intersect.position, lpos_intersect);
914 hitR->detectorElement()->surface(hitR->identify()).globalToLocal(hitR->globalPosition(), hitR->globalPosition(), lPosR);
917 if(std::abs(lpos_intersect.x() - lPosR.x()) > 87.5 )
continue;
923 seeds.emplace_back(std::move(seed));
926 if (usedLayerR) ++seedingLayersR;
929 if (usedLayerL) ++seedingLayersL;
961 NSWSeed& seed,
const std::set<unsigned int>&
exclude,
bool useStereo)
const {
962 ATH_MSG_VERBOSE(
" getClustersOnSegment: layers " << orderedclusters.size());
965 for (
const MeasVec& surfHits : orderedclusters) {
980 ATH_MSG_VERBOSE(
" getClustersOnSegment: returning " << nHitsAdded <<
" hits ");
987 std::vector<NSWSeed> seeds;
989 if (orderedClusters.empty())
return seeds;
991 std::vector<std::vector<const Muon::sTgcPrepData*>> sTgcIP(4);
992 std::vector<std::vector<const Muon::sTgcPrepData*>> sTgcHO(4);
995 for (
int iml : {1, 2}) {
996 int il = (iml == 1) ? 0 : orderedClusters.size() - 1;
997 int iend = (iml == 1) ? orderedClusters.size() : -1;
998 int idir = (iml == 1) ? 1 : -1;
999 unsigned int nLayersWithHitMatch{0};
1002 for (;
il != iend;
il += idir) {
1003 double lastDistance{1000.};
1004 if (nLayersWithHitMatch >= sTgcIP.size()) {
1005 sTgcIP.resize(nLayersWithHitMatch + 1);
1006 sTgcHO.resize(nLayersWithHitMatch + 1);
1008 std::vector<const Muon::sTgcPrepData*>& matchedHits =
1009 (iml == 1) ? sTgcIP.at(nLayersWithHitMatch) : sTgcHO.at(nLayersWithHitMatch);
1014 if (!padHit)
continue;
1020 if (!design)
continue;
1032 double etaDistance = std::abs(padHit->
localPosition().y() - segLocPosOnSurf[1]);
1033 if (etaDistance > 0.5 * chWidth)
continue;
1034 ATH_MSG_DEBUG(
" etaDistance " << etaDistance <<
" between pad center and position on the pad.");
1036 if (matchedHits.empty()) {
1038 matchedHits.push_back(padHit);
1040 }
else if (std::abs(etaDistance - lastDistance) < 0.001) {
1042 matchedHits.push_back(padHit);
1043 ATH_MSG_DEBUG(
" added etaDistance: " << etaDistance <<
" size " << matchedHits.size());
1044 }
else if (etaDistance < lastDistance) {
1046 matchedHits.clear();
1047 matchedHits.push_back(padHit);
1048 ATH_MSG_DEBUG(
" replacing best etaDistance with: " << etaDistance);
1052 lastDistance = etaDistance;
1055 if (!matchedHits.empty()) ++nLayersWithHitMatch;
1060 if (!nLayersWithHitMatch)
return seeds;
1065 std::vector<std::pair<double, double>> sTgcIP_phiRanges =
getPadPhiOverlap(sTgcIP);
1066 std::vector<std::pair<double, double>> sTgcHO_phiRanges =
getPadPhiOverlap(sTgcHO);
1075 for (
const std::pair<double, double>& range1 : sTgcIP_phiRanges) {
1076 double midPhi1 = 0.5 * (range1.first + range1.second);
1079 surfPrdL1.localToGlobal(lp1, gpL1, gpL1);
1087 surfPrdL1.localToGlobal(lp1_lowPhi, gp1_lowPhi, gp1_lowPhi);
1088 surfPrdL1.localToGlobal(lp1_highPhi, gp1_highPhi, gp1_highPhi);
1090 Trk::Intersection intersect_lowPhi = surfPrdL2.straightLineIntersection(gp1_lowPhi, gp1_lowPhi,
false,
false);
1091 Trk::Intersection intersect_highPhi = surfPrdL2.straightLineIntersection(gp1_highPhi, gp1_highPhi,
false,
false);
1096 surfPrdL2.globalToLocal(intersect_lowPhi.
position, intersect_lowPhi.
position, lpIntersect_lowPhi);
1097 surfPrdL2.globalToLocal(intersect_highPhi.
position, intersect_highPhi.
position, lpIntersect_highPhi);
1099 for (
const std::pair<double, double>& range2 : sTgcHO_phiRanges) {
1100 double midPhi2 = 0.5 * (range2.first + range2.second);
1102 if( lpIntersect_highPhi.x() - range2.first < 0 || range2.second - lpIntersect_lowPhi.x() < 0) {
1103 ATH_MSG_DEBUG(
" segmentSeedFromPads: skipping seed with non-overlapping phi ranges");
1109 surfPrdL2.localToGlobal(lp2, gpL2, gpL2);
1113 seeds.emplace_back(
this, gpL1, gDir);
1117 ATH_MSG_DEBUG(
" segmentSeedFromPads: seeds.size() " << seeds.size());
1124 std::vector<NSWSeed> seeds;
1125 std::array<unsigned int, 4>
layers{};
1126 unsigned int trials{0}, used_layers{0};
1128 constexpr
size_t lastMMLay = 11;
1129 std::vector<NSWSeed> laySeeds;
1132 std::stringstream sstr{};
1133 for (
int e4 =
std::min(lastMMLay, orderedClusters.size() -1);
e4 >= 3 ; --
e4) {
1135 for (
int e3 =
e4 -1 ;
e3 >= 2; --
e3) {
1141 const unsigned int old_trials = trials;
1143 if (old_trials == trials)
continue;
1145 used_layers += !laySeeds.empty();
1146 seeds.insert(seeds.end(), std::make_move_iterator(laySeeds.begin()),
1147 std::make_move_iterator(laySeeds.end()));
1150 sstr<<
" Attempts thus far "<<old_trials<<
" attempts now "<<trials<<
" --- "<<
e1<<
","<<
e2<<
","<<
e3<<
","<<
e4<<std::endl;
1152 sstr<<
"Layer: "<<lay<<
" number of measurements "<<orderedClusters[lay].size()<<std::endl;
1154 sstr<<
" **** "<<
print(meas)<<std::endl;
1156 sstr<<std::endl<<std::endl<<std::endl;
1163 if (trials > 100000) {
1166 ATH_MSG_VERBOSE(
"Out of "<<trials<<
" possible seeds, "<<seeds.size()<<
" were finally built. Used in total "<<used_layers<<
" layers");
1170 #if defined(FLATTEN)
1174 std::array<unsigned int,4> selLayers,
1175 unsigned int& trial_counter)
const {
1176 std::vector<NSWSeed> seeds{};
1178 std::array<unsigned int, 4> lay_ord{};
1179 for (
size_t s = 0;
s < selLayers.size(); ++
s) {
1180 unsigned int lay = selLayers[
s];
1188 else lay_ord[
s] = 3;
1190 auto swap_strips = [&selLayers, &lay_ord] (
unsigned int i,
unsigned j){
1195 if (lay_ord[0] == lay_ord[1]){
1196 if (lay_ord[1] != lay_ord[2]) swap_strips(1,2);
1197 else if (lay_ord[1] != lay_ord[3]) swap_strips(1,3);
1204 if (lay_ord[2] == lay_ord[3]) {
1207 if (lay_ord[3] != lay_ord[0] && lay_ord[3] != lay_ord[1]) swap_strips(3,0);
1209 ATH_MSG_VERBOSE(
"No way to rearrange the strips such that the latter two strips cross.");
1214 std::array<SeedMeasurement, 4> base_seed{};
1215 for (
size_t s = 0;
s < selLayers.size(); ++
s) {
1216 unsigned int lay = selLayers[
s];
1217 base_seed[
s] = orderedClusters[lay].front();
1220 const double A = (base_seed[1].pos().z() - base_seed[0].pos().z());
1221 const double G = (base_seed[2].pos().z() - base_seed[0].pos().z());
1222 const double K = (base_seed[3].pos().z() - base_seed[0].pos().z());
1225 diamond.block<2, 1>(0, 1) = ((
A -
G) * base_seed[3].dirDot(base_seed[1]) * base_seed[1].dir() +
G * base_seed[3].dirDot(base_seed[0]) * base_seed[0].dir() -
A * base_seed[3].dir()).block<2, 1>(0, 0);
1226 diamond.block<2, 1>(0, 0) = ((K -
A) * base_seed[2].dirDot(base_seed[1]) * base_seed[1].dir() - K * base_seed[2].dirDot(base_seed[0]) * base_seed[0].dir() +
A * base_seed[2].dir()).block<2, 1>(0, 0);
1228 if (std::abs(diamond.determinant()) < std::numeric_limits<float>::epsilon()) {
1230 << diamond << std::endl
1231 <<
" is singular " << diamond.determinant() <<
" with rank "
1232 << (Eigen::FullPivLU<
AmgSymMatrix(2)>{diamond}.rank()));
1237 << diamond << std::endl
1238 <<
"May give a couple of stereo seeds " << diamond.determinant());
1241 const AmgSymMatrix(2) seed_builder = diamond.inverse();
1243 const double KmG = K-
G;
1244 const double KmA = K-
A;
1245 const double AmG =
A-
G;
1247 const double TwoDotZero = base_seed[2].dirDot(base_seed[0]);
1248 const double ThreeDotZero = base_seed[3].dirDot(base_seed[0]);
1249 const double ThreeDotOne = base_seed[3].dirDot(base_seed[1]);
1250 const double TwoDotOne = base_seed[2].dirDot(base_seed[1]);
1252 auto estimate_muon = [&] () -> std::optional<std::array<double,2>> {
1253 const Amg::Vector3D Y0 = K * base_seed[2].pos() -
G * base_seed[3].pos() - KmG * base_seed[0].pos();
1254 const Amg::Vector3D Y1 = AmG * base_seed[3].pos() - KmG * base_seed[1].pos() + KmA * base_seed[2].pos();
1255 const double Y0dotE0 = base_seed[0].dirDot(Y0);
1256 const double Y1dotE1 = base_seed[1].dirDot(Y1);
1258 const AmgVector(2) centers = (KmG * (base_seed[0].pos() - base_seed[1].pos()) +
1259 Y0dotE0 * base_seed[0].
dir() +
1260 A * (base_seed[3].pos() - base_seed[2].pos()) -
1261 Y1dotE1 * base_seed[1].
dir())
1264 const AmgVector(2) sol_pars = seed_builder * centers;
1265 const std::array<double, 4> lengths{(Y0dotE0 + K * sol_pars[0] * TwoDotZero -
G * sol_pars[1] * ThreeDotZero) / KmG,
1266 (Y1dotE1 + AmG * sol_pars[1] *ThreeDotOne + KmA * sol_pars[0] * TwoDotOne) / KmG, sol_pars[0], sol_pars[1]};
1270 std::optional<Amg::Vector3D> seg_pos{std::nullopt}, seg_dir{std::nullopt};
1271 for (
unsigned int i = 0;
i < base_seed.size(); ++
i) {
1277 seg_pos = std::make_optional<Amg::Vector3D>(base_seed[0].
pos() +
1278 lengths[0] * base_seed[0].
dir());
1282 seg_dir = std::make_optional<Amg::Vector3D>((base_seed[1].
pos() +
1283 lengths[1] *base_seed[1].
dir() - (*seg_pos)).
unit());
1286 std::optional<double> mu_crossing = Amg::intersect<3>(*seg_pos, *seg_dir, base_seed[
i].
pos(),base_seed[
i].
dir());
1288 <<
" ("<< std::string( halfLength > std::abs(lengths[
i]) ?
"inside" :
"outside")<<
" wedge) "
1289 << halfLength <<
" vs. "<<std::abs(lengths[
i])<<
" crossing point: "<<std::abs(*mu_crossing));
1290 }
else if (!
accept)
return std::nullopt;
1292 if (!
accept)
return std::nullopt;
1293 return std::make_optional<std::array<double,2>>({lengths[0], lengths[1]});
1298 MeasVec::const_iterator begin2{orderedClusters[selLayers[1]].begin()};
1299 MeasVec::const_iterator begin3{orderedClusters[selLayers[2]].begin()};
1300 MeasVec::const_iterator begin4{orderedClusters[selLayers[3]].begin()};
1302 const MeasVec::const_iterator end2{orderedClusters[selLayers[1]].end()};
1303 const MeasVec::const_iterator end3{orderedClusters[selLayers[2]].end()};
1304 const MeasVec::const_iterator end4{orderedClusters[selLayers[3]].end()};
1307 base_seed[0] = lay1;
1308 for (MeasVec::const_iterator lay2 = begin2; lay2 != end2; ++lay2) {
1310 base_seed[1] = *lay2;
1323 for (MeasVec::const_iterator lay3 = begin3; lay3 != end3; ++lay3) {
1332 base_seed[2] = *lay3;
1334 for (MeasVec::const_iterator lay4 = begin4 ; lay4 != end4; ++lay4) {
1343 base_seed[3] = (*lay4);
1344 std::optional<std::array<double, 2>> isects = estimate_muon();
1346 if (!isects)
continue;
1347 NSWSeed seed{
this, base_seed, *isects};
1349 if (seed.size() < 4)
continue;
1351 const double eta = std::abs(seed.dir().eta());
1352 if (eta < minEtaNSW || eta > maxEtaNSW) {
1355 if (seed.dir().block<2,1>(0,0).dot(seed.pos().block<2,1>(0,0)) < 0.)
continue;
1358 const double deltaPhi = std::abs(seed.dir().deltaPhi(seed.pos()));
1359 if (
deltaPhi > sector_mapping.sectorWidth(
m_idHelperSvc->sector(base_seed[0]->identify())))
continue;
1361 getClustersOnSegment(orderedClusters, seed, {selLayers[0], selLayers[1],selLayers[2], selLayers[3]});
1362 seeds.emplace_back(std::move(seed));
1381 static const double minTanTheta = 0.75 / std::sinh(maxEtaNSW);
1382 static const double maxTanTheta = 1.25 / std::sinh(minEtaNSW);
1384 const double minDR = minTanTheta * std::abs(dZ);
1385 const double maxDR = maxTanTheta * std::abs(dZ);
1387 const std::pair<double, double> rad1 =
coveredRadii(meas1);
1388 const std::pair<double, double> rad2 =
coveredRadii(meas2);
1390 const double dlR = rad2.first - rad1.first;
1391 const double drR = rad2.second - rad1.second;
1394 <<std::endl<<
". Separation in dR (left/right) "<<dlR<<
"/"<<drR<<
", dZ "<<dZ
1395 <<
" --> dR has to be in "<<minDR<<
" "<<maxDR);
1396 if ((std::abs(dlR) < minDR && std::abs(drR) < minDR) ||
1397 (dZ > 0 && dlR <0 && drR <0) || (dZ < 0 && dlR >0 && drR > 0)) {
1399 }
if (std::abs(dlR) > maxDR && std::abs(drR) > maxDR && dlR * drR > 0.){
1406 const int chNum =
channel(meas);
1412 return std::make_pair(radLeft, radRight);
1415 std::vector<NSWSeed> seeds;
1416 seeds.reserve(unresolved.size());
1417 std::sort(unresolved.begin(), unresolved.end(),[](
const NSWSeed&
a,
const NSWSeed&
b){
1418 return a.chi2() < b.chi2();
1420 for (
NSWSeed& seed : unresolved) {
1421 bool add_seed{
true};
1433 if (add_seed) seeds.push_back(std::move(seed));
1435 ATH_MSG_VERBOSE(seeds.size()<<
" out of "<<unresolved.size()<<
" passed the overlap removal");
1441 const std::vector<std::vector<const Muon::sTgcPrepData*>>& pads)
const {
1446 std::vector<std::vector<double>> padsPhiL, padsPhiR;
1447 std::vector<double> padsPhiC;
1450 for (
const std::vector<const Muon::sTgcPrepData*>& surfHits : pads) {
1452 std::vector<double> surfPadsPhiL, surfPadsPhiR;
1467 bool samePhi = std::find_if(padsPhiC.begin(), padsPhiC.end(), [&hitPadX, &halfWidthX](
const double prevPadPhi) {
1468 return std::abs(hitPadX - prevPadPhi) < 0.9 * halfWidthX;
1469 }) != padsPhiC.end();
1471 if (samePhi)
continue;
1474 surfPadsPhiL.push_back(hitPadX - halfWidthX);
1475 surfPadsPhiR.push_back(hitPadX + halfWidthX);
1476 padsPhiC.push_back(hitPadX);
1480 padsPhiL.push_back(std::move(surfPadsPhiL));
1481 padsPhiR.push_back(std::move(surfPadsPhiR));
1484 unsigned int nSurf = padsPhiR.size();
1488 unsigned int nCombos{1};
1489 for (
const std::vector<double>& surfPadsPhiR : padsPhiR) {
1490 if (!surfPadsPhiR.empty()) nCombos *= surfPadsPhiR.size();
1493 std::vector<std::pair<double, double>> phiOverlap;
1494 phiOverlap.reserve(nCombos);
1496 if (nCombos <= 100) {
1497 unsigned int N{nCombos};
1498 for (
unsigned int isurf{0}; isurf < nSurf; ++isurf) {
1499 if (padsPhiR[isurf].empty())
continue;
1500 unsigned int nSurfHits = padsPhiR[isurf].size();
1503 for (
unsigned int icombo{0}; icombo < nCombos; ++icombo) {
1505 unsigned int padIdx = (icombo /
N) % nSurfHits;
1508 phiOverlap.emplace_back(padsPhiL[isurf][padIdx], padsPhiR[isurf][padIdx]);
1511 phiOverlap[icombo].first =
std::max(padsPhiL[isurf][padIdx], phiOverlap[icombo].
first);
1512 phiOverlap[icombo].second =
std::min(padsPhiR[isurf][padIdx], phiOverlap[icombo].
second);
1518 phiOverlap.erase(std::remove_if(phiOverlap.begin(), phiOverlap.end(),
1519 [](std::pair<double, double>&
range) { return range.first >= range.second; }),
1521 ATH_MSG_DEBUG(
"Pad seeding - #combinations initial: " << nCombos
1522 <<
", after cleaning for non overlapping pads: " << phiOverlap.size());
1526 for (
unsigned int isurf{0}; isurf < nSurf; ++isurf) {
1527 unsigned int nSurfHits = padsPhiR[isurf].size();
1528 for (
unsigned int ihit{0}; ihit < nSurfHits; ++ihit) {
1529 phiOverlap.emplace_back(padsPhiL[isurf][ihit], padsPhiR[isurf][ihit]);
1532 ATH_MSG_DEBUG(
"Pad seeding - #combinations: " << nCombos <<
" is too large. Seeding from" << phiOverlap.size()
1533 <<
" individual pads.");
1547 std::unique_ptr<const Muon::MuonClusterOnTrack> newClus;
1557 calibratedClusters.emplace_back(seed.newCalibClust(std::move(newClus)));
1561 calibratedClusters.emplace_back(seed.newCalibClust(std::move(newClus)));
1565 return calibratedClusters;
1571 std::stringstream sstr{};
1579 std::stringstream sstr{};
1581 <<
" pointing to (" <<
to_string(
cl.dir())<<
" cluster size: "<<clusterSize(
cl);
1586 std::stringstream sstr{};
1588 sstr<<
" *** "<<
print(
cl)<<std::endl;
1593 std::stringstream sstr{};
1594 unsigned int lay{0};
1595 for (
const MeasVec& clusts: sortedVec){
1596 sstr<<
"Clusters in Layer: "<<(lay+1)<<std::endl;
1597 sstr<<
"#################################################"<<std::endl;
1598 sstr<<
print(clusts);
1612 prunedMeas.reserve(clustInLay.size());
1613 const unsigned int firstCh =
channel(clustInLay[0]);
1614 const unsigned int lastCh =
channel(clustInLay[clustInLay.size() -1]);
1616 const unsigned int deltaCh = lastCh - firstCh;
1619 firstCh<<
", highest channel: "<<lastCh<<
" bin width: "<<
m_ocupMmBinWidth<<
" number of bins"<<
nBins);
1622 occupancyHisto.resize(
nBins);
1624 bin.reserve(clustInLay.size());
1630 occupancyHisto[
bin].push_back(std::move(meas));
1640 for (
unsigned int i = 0;
i < occupancyHisto.size() -1; ++
i) {
1642 ATH_MSG_VERBOSE(
"The two neighbouring bins "<<
i<<
"&"<<(
i+1)<<
" have too many clusters "<<std::endl<<
1643 print(occupancyHisto[
i])<<std::endl<<
print(occupancyHisto[
i+1]));
1644 occupancyHisto[
i].clear();
1645 occupancyHisto[
i+1].clear();
1651 std::make_move_iterator(
bin.end()),
1652 std::back_inserter(prunedMeas));
1655 ATH_MSG_VERBOSE(
"Number of measurements before pruning "<<clustInLay.size()<<
" number of measurments survived the pruning "<<prunedMeas.size());