22 #include <GaudiKernel/PhysicalConstants.h>
23 #include <Minuit2/Minuit2Minimizer.h>
24 #include <Math/Minimizer.h>
29 using namespace SegmentFit;
37 copied.reserve(
hits.size());
39 [](
const auto& hit) {
return std::make_unique<MuonR4::CalibratedSpacePoint>(*hit);});
46 toRet.chi2 = toCopy.
chi2;
47 toRet.nDoF = toCopy.
nDoF;
52 std::size_t before =
hits.size();
55 return a->type() == xAOD::UncalibMeasType::Other;
57 return before !=
hits.size();
69 ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
71 m_ambiSolver = std::make_unique<SegmentAmbiSolver>(
name(), std::move(
cfg));
72 return StatusCode::SUCCESS;
77 ATH_CHECK(retrieveContainer(ctx, m_geoCtxKey, gctx));
79 ATH_CHECK(retrieveContainer(ctx, m_seedKey, segmentSeeds));
82 ATH_CHECK(writeSegments.record(std::make_unique<SegmentContainer>()));
83 std::vector<std::unique_ptr<Segment>> allSegments{};
85 std::vector<std::unique_ptr<Segment>> segments = fitSegmentSeed(ctx, *gctx, seed);
86 if (m_visionTool.isEnabled() && segments.size() > 1) {
87 auto drawFinalReco = [
this, &segments, &gctx, &ctx,&seed](
const std::string& nameTag) {
90 segmentLines.push_back(
drawLabel(
std::format(
"# segments: {:d}", segments.size()), 0.12, yLegend, 14));
92 for (
const std::unique_ptr<Segment>& seg : segments) {
97 std::stringstream signStream{};
98 signStream<<
"#chi^{2}/nDoF: "<<
std::format(
"{:.2f}", seg->chi2() / seg->nDoF())<<
"("<<seg->nDoF()<<
"), ";
106 segmentLines.push_back(
drawLabel(signStream.str(), 0.12, yLegend, 12));
110 m_visionTool->visualizeBucket(ctx, *seed->parentBucket(), nameTag, std::move(segmentLines));
112 drawFinalReco(
"all segments");
113 const unsigned int nBeforeAmbi = segments.size();
114 segments = m_ambiSolver->resolveAmbiguity(*gctx, std::move(segments));
115 if (nBeforeAmbi != segments.size()) {
116 drawFinalReco(
"post ambiguity");
119 allSegments.insert(allSegments.end(), std::make_move_iterator(segments.begin()),
120 std::make_move_iterator(segments.end()));
122 resolveAmbiguities(*gctx, allSegments);
123 writeSegments->insert(writeSegments->end(),
124 std::make_move_iterator(allSegments.begin()),
125 std::make_move_iterator(allSegments.end()));
126 ATH_MSG_VERBOSE(
"Found in total "<<writeSegments->size()<<
" segments. ");
127 return StatusCode::SUCCESS;
130 template <
class ContainerType>
133 const ContainerType*& contToPush)
const {
134 contToPush =
nullptr;
137 return StatusCode::SUCCESS;
141 contToPush = readHandle.cptr();
142 return StatusCode::SUCCESS;
150 if (calibHits.empty()) {
156 if (m_doBeamspotConstraint) {
157 unsigned int numPhi =
std::accumulate(calibHits.begin(), calibHits.end(),0,
158 [](
unsigned int n,
const HitVec::value_type& hit){
159 return n + (hit->fitState() == CalibratedSpacePoint::State::Valid &&
163 const Amg::Transform3D globToLoc{calibHits[0]->spacePoint()->msSector()->globalToLocalTrans(gctx)};
167 covariance(0,0) =
std::pow(m_beamSpotR, 2);
168 covariance(1,1) =
std::pow(m_beamSpotR, 2);
169 covariance(2,2) =
std::pow(m_beamSpotL, 2);
171 covariance = jacobian * covariance * jacobian.transpose();
172 AmgSymMatrix(2) beamSpotCov{covariance.block<2,2>(0,0)};
174 beamSpotSP->setCovariance<2>(std::move(beamSpotCov));
176 calibHits.push_back(std::move(beamSpotSP));
180 const Amg::Transform3D& locToGlob{calibHits[0]->spacePoint()->msSector()->localToGlobalTrans(gctx)};
185 fitCfg.doTimeFit = m_doT0Fit;
188 return fitter.fitSegment(ctx, std::move(calibHits), startPars, locToGlob);
191 data.segmentPars = startPars;
194 data.hasPhi = c2f.hasPhiMeas();
195 data.timeFit = c2f.doTimeFit();
196 data.nDoF = c2f.nDoF();
197 if (
data.nDoF <= 0) {
201 ROOT::Minuit2::Minuit2Minimizer minimizer((
name() +
std::to_string(ctx.eventID().event_number())).c_str());
203 minimizer.SetMaxFunctionCalls(100000);
204 minimizer.SetTolerance(0.0001);
205 minimizer.SetPrintLevel(-1);
206 minimizer.SetStrategy(1);
239 minimizer.SetFunction(c2f);
241 if (!minimizer.Minimize() || !minimizer.Hesse()) {
242 data.calibMeasurements = std::move(calibHits);
244 const double* xs = minimizer.X();
245 const double* errs = minimizer.Errors();
248 data.segmentPars[
p] = xs[
p];
249 data.segmentParErrs(
p,
p) = errs[
p];
251 std::optional<double> timeOfArrival{std::nullopt};
252 const auto [locPos,
locDir] =
data.makeLine();
254 timeOfArrival = std::make_optional<double>((locToGlob*locPos).
mag() *
inv_c);
256 data.nIter = minimizer.NCalls();
257 data.calibMeasurements = c2f.release(xs);
260 data.calibMeasurements, msgStream());
261 data.chi2PerMeasurement = std::move(chiPerMeas);
262 data.chi2 = finalChi2;
263 data.converged =
true;
267 std::vector<std::unique_ptr<Segment>>
273 std::vector<std::unique_ptr<Segment>> segments{};
277 genCfg.recalibSeedCircles = m_recalibSeed;
278 genCfg.fastSeedFit = m_refineSeed;
279 genCfg.calibrator = m_calibTool.get();
285 if (m_visionTool.isEnabled()) {
288 while(
auto s = drawMe.nextSeed(ctx)) {
291 seedLines.push_back(
drawLabel(
std::format(
"possible seeds: {:d}", drawMe.numGenerated()), 0.15, 0.85, 14));
292 m_visionTool->visualizeSeed(ctx, *patternSeed,
"pattern", std::move(seedLines));
296 while (
auto seed = seedGen.nextSeed(ctx)) {
298 data.segmentPars = seed->parameters;
299 data.calibMeasurements = std::move(seed->measurements);
301 if (m_visionTool.isEnabled()) {
302 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
303 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Pre fit {:d}", seedGen.numGenerated()));
305 data = fitSegmentHits(ctx, gctx, seed->parameters, std::move(
data.calibMeasurements));
306 data.nIter += seed->nIter;
307 if (m_visionTool.isEnabled() &&
data.converged) {
308 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
309 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Intermediate fit {:d}", seedGen.numGenerated()));
311 if (!removeOutliers(ctx, gctx, *patternSeed,
data)) {
314 if (!plugHoles(ctx, gctx, *patternSeed,
data)) {
317 if (m_visionTool.isEnabled()) {
318 auto seedCopy = convertToSegment(locToGlob, patternSeed,
copy(
data));
319 m_visionTool->visualizeSegment(ctx, *seedCopy,
std::format(
"Final fit {:d}", seedGen.numGenerated()));
321 segments.push_back(convertToSegment(locToGlob, patternSeed, std::move(
data)));
328 const auto [locPos,
locDir] =
data.makeLine();
333 auto finalSeg = std::make_unique<Segment>(std::move(globPos), std::move(globDir),
334 patternSeed, std::move(
data.calibMeasurements),
336 finalSeg->setCallsToConverge(
data.nIter);
337 finalSeg->setChi2PerMeasurement(std::move(
data.chi2PerMeasurement));
338 finalSeg->setParUncertainties(std::move(
data.segmentParErrs));
354 if (
data.nDoF<=0 ||
data.calibMeasurements.empty()) {
355 ATH_MSG_VERBOSE(
"No degree of freedom available. What shall be removed?!");
359 const auto [segPos, segDir] =
data.makeLine();
361 if (
data.converged &&
data.chi2 /
data.nDoF < m_outlierRemovalCut) {
364 <<
". Don't remove outliers");
376 std::sort(
data.calibMeasurements.begin(),
data.calibMeasurements.end(),
377 [&,
this](
const HitVec::value_type&
a,
const HitVec::value_type&
b){
378 return SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *a, msgStream()) <
379 SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *b, msgStream());
384 data.nDoF -=
data.calibMeasurements.back()->measuresEta();
385 data.nDoF -=
data.calibMeasurements.back()->measuresPhi();
386 if (m_doT0Fit &&
data.calibMeasurements.back()->measuresTime() &&
391 std::vector<HoughHitType> uncalib{};
392 for (
const HitVec::value_type&
calib :
data.calibMeasurements) {
393 uncalib.push_back(
calib->spacePoint());
398 data = std::move(newAttempt);
399 if (m_visionTool.isEnabled()) {
400 auto seedCopy = convertToSegment(seed.msSector()->localToGlobalTrans(gctx), &seed,
copy(
data));
401 m_visionTool->visualizeSegment(ctx, *seedCopy,
"Bad fit recovery");
406 return removeOutliers(ctx, gctx, seed,
data);
415 <<
", chi2: "<<beforeRecov.
chi2<<
", nDoF: "<<beforeRecov.
nDoF<<
", "
418 std::unordered_set<const SpacePoint*> usedSpacePoint{};
420 usedSpacePoint.insert(hit->spacePoint());
425 bool hasCandidate{
false};
427 for (
const std::vector<HoughHitType>& mdtLayer : hitLayers.mdtHits()) {
430 if (usedSpacePoint.count(mdtHit)) {
433 const double dist =
Amg::lineDistance(locPos,
locDir, mdtHit->positionInChamber(), mdtHit->directionInChamber());
435 if (dist >= dc->readoutElement()->innerTubeRadius()) {
440 ATH_MSG_VERBOSE(__func__<<
"() :"<<__LINE__<<
" Candidate hit for recovery "<<m_idHelperSvc->toString(mdtHit->identify())<<
", chi2: "<<
pull);
441 if (
pull <= m_recoveryPull) {
443 candidateHits.push_back(std::move(calibHit));
446 candidateHits.push_back(std::move(calibHit));
453 std::make_move_iterator(candidateHits.begin()),
454 std::make_move_iterator(candidateHits.end()));
455 eraseWrongHits(gctx, beforeRecov);
456 return beforeRecov.
nDoF > 0;
461 if (m_doBeamspotConstraint) {
465 candidateHits.insert(candidateHits.end(), std::make_move_iterator(copied.begin()), std::make_move_iterator(copied.end()));
473 for (HitVec::value_type& hit : copiedCandidates) {
477 eraseWrongHits(gctx, beforeRecov);
483 if (redChi2 < m_outlierRemovalCut || (beforeRecov.
nDoF == 0) || redChi2 < beforeRecov.
chi2 / beforeRecov.
nDoF) {
486 beforeRecov = std::move(recovered);
489 bool runAnotherTrial =
false;
491 if (m_doBeamspotConstraint) {
498 runAnotherTrial =
true;
501 if (!runAnotherTrial) {
504 recovered = fitSegmentHits(ctx, gctx, beforeRecov.
segmentPars, std::move(copied));
512 if (redChi2 < m_outlierRemovalCut || redChi2 < beforeRecov.
chi2 / beforeRecov.
nDoF) {
514 beforeRecov = std::move(recovered);
520 eraseWrongHits(gctx, beforeRecov);
522 for (HitVec::value_type& hit : copiedCandidates) {
532 [&segPos, &segDir](
const HitVec::value_type& hit){
533 if (hit->fitState() != CalibratedSpacePoint::State::Outlier) {
538 const double dist = Amg::lineDistance(segPos, segDir, hit->positionInChamber(), hit->directionInChamber());
539 const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(hit->spacePoint()->primaryMeasurement());
540 return dist >= dc->readoutElement()->innerTubeRadius();
543 }), candidate.calibMeasurements.end());
545 for (
const auto& hit : candidate.calibMeasurements) {
547 chamber = hit->spacePoint()->msSector();
551 std::optional<double> timeOfFlight = candidate.timeFit ? std::make_optional<double>((
chamber->localToGlobalTrans(gctx)*segPos).mag() *
inv_c)
554 candidate.calibMeasurements, msgStream());
555 ATH_MSG_VERBOSE(
"The measurements before "<<candidate.chi2PerMeasurement<<
", after: "<<updatedMeasChi2<<
", chi2: "<<updateChi2);
556 candidate.chi2PerMeasurement = std::move(updatedMeasChi2);
559 std::vector<std::unique_ptr<Segment>>& segmentCandidates)
const {
560 using SegmentVec = std::vector<std::unique_ptr<Segment>>;
561 ATH_MSG_VERBOSE(
"Resolve ambiguities amongst "<<segmentCandidates.size()<<
" segment candidates. ");
562 std::unordered_map<const MuonGMR4::SpectrometerSector*, SegmentVec> candidatesPerChamber{};
564 for (std::unique_ptr<Segment>& sortMe : segmentCandidates) {
566 candidatesPerChamber[chamb].push_back(std::move(sortMe));
568 segmentCandidates.clear();
569 for (
auto& [
chamber, resolveMe] : candidatesPerChamber) {
570 SegmentVec resolvedSegments = m_ambiSolver->resolveAmbiguity(gctx, std::move(resolveMe));
571 segmentCandidates.insert(segmentCandidates.end(),
572 std::make_move_iterator(resolvedSegments.begin()),
573 std::make_move_iterator(resolvedSegments.end()));