15 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
16 #include "Acts/MagneticField/MagneticFieldContext.hpp"
17 #include "Acts/Surfaces/PerigeeSurface.hpp"
18 #include "Acts/Utilities/CalibrationContext.hpp"
19 #include "Acts/Definitions/TrackParametrization.hpp"
20 #include "Acts/Definitions/Units.hpp"
21 #include "Acts/Propagator/EigenStepper.hpp"
22 #include "Acts/Surfaces/Surface.hpp"
23 #include "Acts/TrackFitting/GsfMixtureReduction.hpp"
24 #include "Acts/Utilities/Helpers.hpp"
25 #include "Acts/EventData/VectorMultiTrajectory.hpp"
26 #include "Acts/EventData/VectorTrackContainer.hpp"
27 #include "Acts/EventData/TrackParameters.hpp"
28 #include "Acts/Utilities/Logger.hpp"
29 #include "Acts/Utilities/CalibrationContext.hpp"
39 #include "Acts/Propagator/DirectNavigator.hpp"
45 #include <type_traits>
46 #include <system_error>
52 const IInterface*
p) :
70 auto field = std::make_shared<ATLASMagneticFieldWrapper>();
71 Acts::MultiEigenStepperLoop<> stepper(
field);
73 logger().cloneWithSuffix(
"Navigator") );
74 Acts::Propagator<Acts::MultiEigenStepperLoop<>,
Acts::Navigator> propagator(std::move(stepper),
76 logger().cloneWithSuffix(
"Prop"));
78 Acts::AtlasBetheHeitlerApprox<6, 5> bha = Acts::makeDefaultBetheHeitlerApprox();
79 m_fitter = std::make_unique<Fitter>(std::move(propagator), std::move(bha),
80 logger().cloneWithSuffix(
"GaussianSumFitter"));
84 Acts::DirectNavigator directNavigator(
logger().cloneWithSuffix(
"DirectNavigator") );
85 Acts::MultiEigenStepperLoop<> stepperDirect(
field);
86 Acts::Propagator<Acts::MultiEigenStepperLoop<>, Acts::DirectNavigator> directPropagator(std::move(stepperDirect),
87 std::move(directNavigator),
88 logger().cloneWithSuffix(
"DirectPropagator"));
89 Acts::AtlasBetheHeitlerApprox<6, 5> bhaDirect = Acts::makeDefaultBetheHeitlerApprox();
90 m_directFitter = std::make_unique<DirectFitter>(std::move(directPropagator),std::move(bhaDirect),
91 logger().cloneWithSuffix(
"DirectGaussianSumFitter"));
95 m_gsfExtensions.updater.connect<&ActsTrk::detail::FitterHelperFunctions::gainMatrixUpdate<ActsTrk::MutableTrackStateBackend>>();
96 m_calibrator = std::make_unique<ActsTrk::detail::TrkMeasurementCalibrator>();
97 m_gsfExtensions.calibrator.connect<&ActsTrk::detail::TrkMeasurementCalibrator::calibrate<ActsTrk::MutableTrackStateBackend>>(
m_calibrator.get());
101 m_gsfExtensions.mixtureReducer.connect<&Acts::reduceMixtureWithKLDistance>();
119 return StatusCode::SUCCESS;
124 std::unique_ptr<Trk::Track>
132 std::unique_ptr<Trk::Track>
track =
nullptr;
133 ATH_MSG_VERBOSE (
"--> enter GaussianSumFitter::fit(Track,,) with Track from author = "
138 ATH_MSG_DEBUG(
"called to refit empty track or track with too little information, reject fit");
144 ATH_MSG_DEBUG(
"input fails to provide track parameters for seeding the GSF, reject fit");
149 std::shared_ptr<Acts::PerigeeSurface> pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
153 const Acts::MagneticFieldContext mfContext =
m_extrapolationTool->getMagneticFieldContext(ctx);
157 Acts::GsfOptions<ActsTrk::MutableTrackStateBackend>
167 std::vector<Acts::SourceLink> trackSourceLinks =
m_ATLASConverterTool->trkTrackToSourceLinks(inputTrack);
179 std::unique_ptr<Trk::Track>
186 std::unique_ptr<Trk::Track>
track =
nullptr;
188 if (inputMeasSet.size() < 2) {
189 ATH_MSG_DEBUG(
"called to refit empty measurement set or a measurement set with too little information, reject fit");
194 std::shared_ptr<Acts::PerigeeSurface> pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
198 const Acts::MagneticFieldContext mfContext =
m_extrapolationTool->getMagneticFieldContext(ctx);
202 Acts::GsfOptions<ActsTrk::MutableTrackStateBackend>
209 gsfOptions.abortOnError =
false;
211 std::vector< Acts::SourceLink > trackSourceLinks;
214 const auto initialParams =
m_ATLASConverterTool->trkTrackParametersToActsParameters(estimatedStartParameters, tgContext);
218 std::vector<const Acts::Surface*> surfaces;
219 surfaces.reserve(inputMeasSet.size());
239 std::unique_ptr<Trk::Track>
254 std::unique_ptr<Trk::Track>
261 ATH_MSG_VERBOSE (
"--> enter GaussianSumFitter::fit(Track,Meas'BaseSet,,)");
265 if (addMeasColl.empty()) {
266 ATH_MSG_DEBUG(
"client tries to add an empty MeasurementSet to the track fit." );
267 return fit(ctx,inputTrack);
272 ATH_MSG_DEBUG(
"called to refit empty track or track with too little information, reject fit");
278 ATH_MSG_DEBUG(
"input fails to provide track parameters for seeding the GSF, reject fit");
282 std::unique_ptr<Trk::Track>
track =
nullptr;
285 std::shared_ptr<Acts::PerigeeSurface> pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
Acts::Vector3::Zero());
288 const Acts::MagneticFieldContext mfContext =
m_extrapolationTool->getMagneticFieldContext(ctx);
292 Acts::GsfOptions<ActsTrk::MutableTrackStateBackend>
298 std::vector<Acts::SourceLink> trackSourceLinks =
m_ATLASConverterTool->trkTrackToSourceLinks(inputTrack);
310 std::unique_ptr<Trk::Track>
318 ATH_MSG_DEBUG(
"Fit of Track with additional PrepRawDataSet not yet implemented");
324 std::unique_ptr<Trk::Track>
337 ATH_MSG_DEBUG(
"input #2 is empty try to fit track 1 alone" );
338 return fit(ctx,intrk1);
343 ATH_MSG_DEBUG(
"input #1 is empty try to fit track 2 alone" );
344 return fit(ctx,intrk2);
349 ATH_MSG_DEBUG(
"input #1 fails to provide track parameters for seeding the GSF, reject fit");
353 std::unique_ptr<Trk::Track>
track =
nullptr;
356 std::shared_ptr<Acts::PerigeeSurface> pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
360 const Acts::MagneticFieldContext mfContext =
m_extrapolationTool->getMagneticFieldContext(ctx);
364 Acts::GsfOptions<ActsTrk::MutableTrackStateBackend>
370 std::vector<Acts::SourceLink> trackSourceLinks =
m_ATLASConverterTool->trkTrackToSourceLinks(intrk1);
371 std::vector<Acts::SourceLink> trackSourceLinks2 =
m_ATLASConverterTool->trkTrackToSourceLinks(intrk2);
372 trackSourceLinks.insert(trackSourceLinks.end(), trackSourceLinks2.begin(), trackSourceLinks2.end());
383 std::unique_ptr< ActsTrk::MutableTrackContainer >
385 const Acts::BoundTrackParameters& ,
386 const Acts::GeometryContext& ,
387 const Acts::MagneticFieldContext& ,
388 const Acts::CalibrationContext& ,
389 const Acts::Surface& )
const
391 ATH_MSG_VERBOSE(
"ACTS seed refit is not implemented in GaussianSumFitterTool");
395 std::unique_ptr< ActsTrk::MutableTrackContainer >
397 const Acts::BoundTrackParameters& ,
398 const Acts::GeometryContext& ,
399 const Acts::MagneticFieldContext& ,
400 const Acts::CalibrationContext& ,
401 const Acts::Surface* )
const
403 ATH_MSG_VERBOSE(
"ACTS uncalib slink refit is not implemented in GaussianSumFitterTool");
409 const EventContext& ctx,
410 const ActsTrk::TrackContainer::ConstTrackProxy&
track,
412 const Acts::PerigeeSurface& pSurface)
const {
415 const Acts::BoundTrackParameters initialParams =
track.createParametersAtReference();
416 std::vector<Acts::SourceLink> sourceLinks;
417 std::vector<const Acts::Surface*> surfSequence;
419 for (
auto ts :
track.trackStates()){
420 surfSequence.push_back(&
ts.referenceSurface());
421 if (!
ts.hasCalibrated()) {
425 if (
ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
430 if (sourceLinks.size() < 2) {
431 ATH_MSG_DEBUG(
"called to refit 0 or 1 sourceLink with too little information, reject fit");
432 return StatusCode::SUCCESS;
437 Acts::CalibrationContext calContext{};
439 Acts::GsfOptions<ActsTrk::MutableTrackStateBackend> gsfOptions =
prepareOptions(tgContext, mfContext, calContext, pSurface);
442 if (!actsTrackingGeometry) {
444 return StatusCode::FAILURE;
452 gsfOptions.extensions = gsfExtensions;
453 gsfOptions.abortOnError =
false;
461 return StatusCode::SUCCESS;
464 std::unique_ptr<Trk::Track>
466 const Acts::GeometryContext& tgContext,
468 Acts::Result<typename ActsTrk::MutableTrackContainer::TrackProxy, std::error_code>& fitResult)
const
470 if (not fitResult.ok())
473 std::unique_ptr<Trk::Track> newtrack =
nullptr;
475 const auto& acts_track = fitResult.value();
476 auto finalTrajectory = std::make_unique<Trk::TrackStates>();
478 int numberOfDeadPixel = 0;
479 int numberOfDeadSCT = 0;
481 std::vector<std::unique_ptr<const Acts::BoundTrackParameters>> actsSmoothedParam;
483 tracks.trackStateContainer().visitBackwards(acts_track.tipIndex(),
484 [&] (
const auto &state) ->
void
487 auto flag = state.typeFlags();
488 const auto* associatedDetEl = state.referenceSurface().associatedDetectorElement();
489 if (not associatedDetEl)
492 const auto* actsElement = dynamic_cast<const ActsDetectorElement*>(associatedDetEl);
496 const auto* upstreamDetEl = actsElement->upstreamDetectorElement();
497 if (not upstreamDetEl)
500 ATH_MSG_VERBOSE(
"Try casting to TRT for if");
501 if (dynamic_cast<const InDetDD::TRT_BaseElement*>(upstreamDetEl))
504 const auto* trkDetElem = dynamic_cast<const Trk::TrkDetElementBase*>(upstreamDetEl);
508 ATH_MSG_VERBOSE(
"trkDetElem type: " << static_cast<std::underlying_type_t<Trk::DetectorElemType>>(trkDetElem->detectorType()));
510 ATH_MSG_VERBOSE(
"Try casting to SiDetectorElement");
511 const auto* detElem = dynamic_cast<const InDetDD::SiDetectorElement*>(upstreamDetEl);
514 ATH_MSG_VERBOSE(
"detElem = " << detElem);
517 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
518 std::unique_ptr<Trk::TrackParameters> parm;
521 if (flag.test(Acts::TrackStateFlag::HoleFlag)){
522 ATH_MSG_VERBOSE(
"State is a hole (no associated measurement), use predicted parameters");
523 const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
525 state.predictedCovariance(),
526 acts_track.particleHypothesis());
527 parm = m_ATLASConverterTool->actsTrackParametersToTrkParameters(actsParam, tgContext);
528 auto boundaryCheck = m_boundaryCheckTool->boundaryCheck(*parm);
530 ATH_MSG_VERBOSE(
"Check if this is a hole, a dead sensors or a state outside the sensor boundary");
531 if(boundaryCheck == Trk::BoundaryCheckResult::DeadElement){
532 if (detElem->isPixel()) {
535 else if (detElem->isSCT()) {
540 } else if (boundaryCheck != Trk::BoundaryCheckResult::Candidate){
544 typePattern.set(Trk::TrackStateOnSurface::Hole);
547 else if (
flag.test(Acts::TrackStateFlag::OutlierFlag) or not state.hasSmoothed()) {
548 ATH_MSG_VERBOSE(
"The state was tagged as an outlier or was missed in the reverse filtering, use filtered parameters");
549 const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
551 state.filteredCovariance(),
552 acts_track.particleHypothesis());
553 parm = m_ATLASConverterTool->actsTrackParametersToTrkParameters(actsParam, tgContext);
554 typePattern.set(Trk::TrackStateOnSurface::Outlier);
558 ATH_MSG_VERBOSE(
"The state is a measurement state, use smoothed parameters");
560 const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
562 state.smoothedCovariance(),
563 acts_track.particleHypothesis());
565 actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam)));
566 parm = m_ATLASConverterTool->actsTrackParametersToTrkParameters(actsParam, tgContext);
567 typePattern.set(Trk::TrackStateOnSurface::Measurement);
570 std::unique_ptr<Trk::MeasurementBase> measState;
571 if (state.hasUncalibratedSourceLink()){
572 auto sl = state.getUncalibratedSourceLink().template get<ATLASSourceLink>();
574 measState = sl->uniqueClone();
576 double nDoF = state.calibratedSize();
581 ATH_MSG_VERBOSE(
"State succesfully creates, adding it to the trajectory");
582 finalTrajectory->insert(finalTrajectory->begin(), perState);
587 const Acts::BoundTrackParameters actsPer(acts_track.referenceSurface().getSharedPtr(),
588 acts_track.parameters(),
589 acts_track.covariance(),
590 acts_track.particleHypothesis());
591 std::unique_ptr<Trk::TrackParameters> per = m_ATLASConverterTool->actsTrackParametersToTrkParameters(actsPer, tgContext);
592 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
595 if (perState) finalTrajectory->insert(finalTrajectory->begin(), perState);
600 newtrack = std::make_unique<Trk::Track>(newInfo, std::move(finalTrajectory),
nullptr);
601 if (newtrack && !m_refitOnly) {
603 if (!newtrack->trackSummary()) {
604 newtrack->setTrackSummary(std::make_unique<Trk::TrackSummary>());
611 m_trkSummaryTool->updateTrackSummary(ctx, *newtrack,
true);
617 const Acts::GsfExtensions<typename ActsTrk::MutableTrackStateBackend>&
618 GaussianSumFitterTool::getExtensions()
const
620 return m_gsfExtensions;
630 Acts::GsfOptions<typename ActsTrk::MutableTrackStateBackend>
631 GaussianSumFitterTool::prepareOptions(
const Acts::GeometryContext& tgContext,
632 const Acts::MagneticFieldContext& mfContext,
633 const Acts::CalibrationContext& calContext,
634 const Acts::PerigeeSurface& surface)
const
636 Acts::PropagatorPlainOptions propagationOption(tgContext, mfContext);
637 propagationOption.maxSteps = m_option_maxPropagationStep;
639 Acts::GsfOptions<typename ActsTrk::MutableTrackStateBackend> gsfOptions(tgContext, mfContext, calContext);
640 gsfOptions.extensions=m_gsfExtensions;
641 gsfOptions.propagatorPlainOptions=propagationOption;
642 gsfOptions.referenceSurface = &surface;
645 gsfOptions.abortOnError =
false;
646 gsfOptions.maxComponents = m_maxComponents;
647 gsfOptions.weightCutoff = m_weightCutOff;
648 gsfOptions.componentMergeMethod = m_componentMergeMethod;