22 #include <GaudiKernel/PhysicalConstants.h>
23 #include <Minuit2/Minuit2Minimizer.h>
24 #include <Math/Minimizer.h>
31 using namespace SegmentFit;
36 copied.reserve(
hits.size());
38 [](
const auto& hit) {
return std::make_unique<MuonR4::CalibratedSpacePoint>(*hit);});
45 toRet.chi2 = toCopy.
chi2;
46 toRet.nDoF = toCopy.
nDoF;
51 std::size_t before =
hits.size();
54 return a->type() == xAOD::UncalibMeasType::Other;
56 return before !=
hits.size();
68 ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
70 m_ambiSolver = std::make_unique<SegmentAmbiSolver>(
name(), std::move(
cfg));
71 return StatusCode::SUCCESS;
80 ATH_CHECK(writeSegments.record(std::make_unique<SegmentContainer>()));
82 std::vector<std::unique_ptr<Segment>> allSegments{};
84 std::vector<std::unique_ptr<Segment>> segments = fitSegmentSeed(ctx, *gctx, seed);
85 if (m_visionTool.isEnabled() && segments.size() > 1) {
86 auto drawFinalReco = [
this, &segments, &gctx, &ctx,&seed](
const std::string& nameTag) {
89 segmentLines.push_back(
drawLabel(
std::format(
"# segments: {:d}", segments.size()), 0.2, yLegend, 14));
91 for (
const std::unique_ptr<Segment>& seg : segments) {
96 std::stringstream signStream{};
97 signStream<<
std::format(
"#chi^{{2}}/nDoF: {:.2f} ({:}), ", seg->chi2() / seg->nDoF(), seg->nDoF());
105 segmentLines.push_back(
drawLabel(signStream.str(), 0.2, yLegend, 13));
109 m_visionTool->visualizeBucket(ctx, *seed->parentBucket(), nameTag, std::move(segmentLines));
111 drawFinalReco(
"all segments");
112 const unsigned int nBeforeAmbi = segments.size();
113 segments = m_ambiSolver->resolveAmbiguity(*gctx, std::move(segments));
114 if (nBeforeAmbi != segments.size()) {
115 drawFinalReco(
"post ambiguity");
117 }
else if (m_visionTool.isEnabled() && segments.empty() &&
118 std::ranges::count_if(seed->getHitsInMax(),[
this](
const SpacePoint* hit){
119 return m_visionTool->isLabeled(*hit);
121 m_visionTool->visualizeSeed(ctx, *seed,
"Failed fit");
123 allSegments.insert(allSegments.end(), std::make_move_iterator(segments.begin()),
124 std::make_move_iterator(segments.end()));
126 resolveAmbiguities(*gctx, allSegments);
127 writeSegments->insert(writeSegments->end(),
128 std::make_move_iterator(allSegments.begin()),
129 std::make_move_iterator(allSegments.end()));
130 ATH_MSG_VERBOSE(
"Found in total "<<writeSegments->size()<<
" segments. ");
131 return StatusCode::SUCCESS;
139 if (calibHits.empty()) {
145 if (m_doBeamspotConstraint) {
146 unsigned int numPhi =
std::accumulate(calibHits.begin(), calibHits.end(),0,
147 [](
unsigned int n,
const HitVec::value_type& hit){
148 return n + (hit->fitState() == CalibratedSpacePoint::State::Valid &&
152 const Amg::Transform3D globToLoc{calibHits[0]->spacePoint()->msSector()->globalToLocalTrans(gctx)};
156 covariance(0,0) =
std::pow(m_beamSpotR, 2);
157 covariance(1,1) =
std::pow(m_beamSpotR, 2);
158 covariance(2,2) =
std::pow(m_beamSpotL, 2);
160 covariance = jacobian * covariance * jacobian.transpose();
161 AmgSymMatrix(2) beamSpotCov{covariance.block<2,2>(0,0)};
163 beamSpotSP->setCovariance<2>(std::move(beamSpotCov));
165 calibHits.push_back(std::move(beamSpotSP));
169 const Amg::Transform3D& locToGlob{calibHits[0]->spacePoint()->msSector()->localToGlobalTrans(gctx)};
173 fitCfg.doTimeFit = m_doT0Fit;
174 fitCfg.reCalibrate = m_recalibInFit;
177 return fitter.fitSegment(ctx, std::move(calibHits), startPars, locToGlob);
179 std::vector<std::unique_ptr<Segment>>
185 std::vector<std::unique_ptr<Segment>> segments{};
189 genCfg.recalibSeedCircles = m_recalibSeed;
190 genCfg.fastSeedFit = m_refineSeed;
191 genCfg.fastSegFitWithT0 = m_doT0Fit;
192 genCfg.calibrator = m_calibTool.get();
198 if (m_visionTool.isEnabled()) {
201 while(
auto s = drawMe.nextSeed(ctx)) {
204 seedLines.push_back(
drawLabel(
std::format(
"possible seeds: {:d}", drawMe.numGenerated()), 0.2, 0.85, 14));
205 m_visionTool->visualizeSeed(ctx, *patternSeed,
"pattern", std::move(seedLines));
209 while (
auto seed = seedGen.nextSeed(ctx)) {
211 data.segmentPars = seed->parameters;
212 data.calibMeasurements = std::move(seed->measurements);
214 if (m_visionTool.isEnabled()) {
215 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
216 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Pre fit {:d}", seedGen.numGenerated()));
218 data = fitSegmentHits(ctx, gctx, seed->parameters, std::move(
data.calibMeasurements));
219 data.nIter += seed->nIter;
220 if (m_visionTool.isEnabled() &&
data.converged) {
221 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
222 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Intermediate fit {:d}", seedGen.numGenerated()));
224 if (!removeOutliers(ctx, gctx, *patternSeed,
data)) {
227 if (!plugHoles(ctx, gctx, *patternSeed,
data)) {
230 if (m_visionTool.isEnabled()) {
231 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
232 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Final fit {:d}", seedGen.numGenerated()));
234 segments.push_back(convertToSegment(locToGlob, patternSeed, std::move(
data)));
241 const auto [locPos,
locDir] =
data.makeLine();
246 auto finalSeg = std::make_unique<Segment>(std::move(globPos), std::move(globDir),
247 patternSeed, std::move(
data.calibMeasurements),
249 finalSeg->setCallsToConverge(
data.nIter);
250 finalSeg->setParUncertainties(std::move(
data.segmentParErrs));
263 if (
data.nDoF<=0 ||
data.calibMeasurements.empty() ||
data.nPrecMeas < m_precHitCut) {
264 ATH_MSG_VERBOSE(
"No degree of freedom available. What shall be removed?!. nDoF: "
265 <<
data.nDoF<<
", n-meas: "<<
data.calibMeasurements);
269 const auto [segPos, segDir] =
data.makeLine();
271 if (
data.converged &&
data.chi2 /
data.nDoF < m_outlierRemovalCut) {
274 <<
". Don't remove outliers");
286 std::ranges::sort(
data.calibMeasurements,
287 [&,
this](
const HitVec::value_type&
a,
const HitVec::value_type&
b){
288 return SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *a, msgStream()) <
289 SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *b, msgStream());
294 data.nDoF -=
data.calibMeasurements.back()->measuresEta();
295 data.nDoF -=
data.calibMeasurements.back()->measuresPhi();
296 if (m_doT0Fit &&
data.calibMeasurements.back()->measuresTime() &&
301 std::vector<HoughHitType> uncalib{};
302 for (
const HitVec::value_type&
calib :
data.calibMeasurements) {
303 uncalib.push_back(
calib->spacePoint());
308 data = std::move(newAttempt);
309 if (m_visionTool.isEnabled()) {
310 auto seedCopy = convertToSegment(seed.msSector()->localToGlobalTrans(gctx), &seed,
copy(
data));
311 m_visionTool->visualizeSegment(ctx, *seedCopy,
"Bad fit recovery");
316 return removeOutliers(ctx, gctx, seed,
data);
325 <<
", chi2: "<<beforeRecov.
chi2<<
", nDoF: "<<beforeRecov.
nDoF<<
", "
328 std::unordered_set<const SpacePoint*> usedSpacePoint{};
330 usedSpacePoint.insert(hit->spacePoint());
335 bool hasCandidate{
false};
337 for (
const std::vector<HoughHitType>& mdtLayer : hitLayers.mdtHits()) {
340 if (usedSpacePoint.count(mdtHit)) {
343 const double dist =
Amg::lineDistance(locPos,
locDir, mdtHit->positionInChamber(), mdtHit->directionInChamber());
345 if (dist >= dc->readoutElement()->innerTubeRadius()) {
350 ATH_MSG_VERBOSE(__func__<<
"() :"<<__LINE__<<
" Candidate hit for recovery "<<m_idHelperSvc->toString(mdtHit->identify())<<
", chi2: "<<
pull);
351 if (
pull <= m_recoveryPull) {
353 candidateHits.push_back(std::move(calibHit));
356 candidateHits.push_back(std::move(calibHit));
363 std::make_move_iterator(candidateHits.begin()),
364 std::make_move_iterator(candidateHits.end()));
365 eraseWrongHits(beforeRecov);
366 return beforeRecov.
nDoF > 0;
371 if (m_doBeamspotConstraint) {
375 candidateHits.insert(candidateHits.end(), std::make_move_iterator(copied.begin()), std::make_move_iterator(copied.end()));
383 for (HitVec::value_type& hit : copiedCandidates) {
387 eraseWrongHits(beforeRecov);
393 if (redChi2 < m_outlierRemovalCut || (beforeRecov.
nDoF == 0) || redChi2 < beforeRecov.
chi2 / beforeRecov.
nDoF) {
396 beforeRecov = std::move(recovered);
399 bool runAnotherTrial =
false;
401 if (m_doBeamspotConstraint) {
404 const auto [beforePos, beforeDir] = beforeRecov.
makeLine();
405 for (HitVec::value_type& copyHit : copied) {
412 std::nullopt, *copyHit, msgStream())) < m_recoveryPull;
414 runAnotherTrial =
true;
419 if (!runAnotherTrial) {
422 recovered = fitSegmentHits(ctx, gctx, beforeRecov.
segmentPars, std::move(copied));
430 if (redChi2 < m_outlierRemovalCut || redChi2 < beforeRecov.
chi2 / beforeRecov.
nDoF) {
432 beforeRecov = std::move(recovered);
438 eraseWrongHits(beforeRecov);
440 for (HitVec::value_type& hit : copiedCandidates) {
450 [&segPos, &segDir](
const HitVec::value_type& hit){
451 if (hit->fitState() != CalibratedSpacePoint::State::Outlier) {
456 const double dist = Amg::lineDistance(segPos, segDir, hit->positionInChamber(), hit->directionInChamber());
457 const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(hit->spacePoint()->primaryMeasurement());
458 return dist >= dc->readoutElement()->innerTubeRadius();
461 }), candidate.calibMeasurements.end());
464 return a->positionInChamber().z() < b->positionInChamber().z();
468 std::vector<std::unique_ptr<Segment>>& segmentCandidates)
const {
469 using SegmentVec = std::vector<std::unique_ptr<Segment>>;
470 ATH_MSG_VERBOSE(
"Resolve ambiguities amongst "<<segmentCandidates.size()<<
" segment candidates. ");
471 std::unordered_map<const MuonGMR4::SpectrometerSector*, SegmentVec> candidatesPerChamber{};
473 for (std::unique_ptr<Segment>& sortMe : segmentCandidates) {
475 candidatesPerChamber[chamb].push_back(std::move(sortMe));
477 segmentCandidates.clear();
478 for (
auto& [
chamber, resolveMe] : candidatesPerChamber) {
479 SegmentVec resolvedSegments = m_ambiSolver->resolveAmbiguity(gctx, std::move(resolveMe));
480 segmentCandidates.insert(segmentCandidates.end(),
481 std::make_move_iterator(resolvedSegments.begin()),
482 std::make_move_iterator(resolvedSegments.end()));