diff --git a/InnerDetector/InDetExample/InDetRecExample/python/TrackingCommon.py b/InnerDetector/InDetExample/InDetRecExample/python/TrackingCommon.py
index 433a8c7a9ed6628e9f9c126cab646262aaf6abe2..99383e9fc4691910df6388e9598cc7969e0a9354 100644
--- a/InnerDetector/InDetExample/InDetRecExample/python/TrackingCommon.py
+++ b/InnerDetector/InDetExample/InDetRecExample/python/TrackingCommon.py
@@ -654,8 +654,8 @@ def getInDetGsfMaterialUpdator(name='InDetGsfMaterialUpdator', **kwargs) :
if 'MaximumNumberOfComponents' not in kwargs :
kwargs=setDefaults(kwargs, MaximumNumberOfComponents = 12)
- from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__GsfMaterialMixtureConvolution
- return Trk__GsfMaterialMixtureConvolution (name = the_name, **kwargs)
+ from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__GsfMaterialMixtureConvolutionLM
+ return Trk__GsfMaterialMixtureConvolutionLM (name = the_name, **kwargs)
@makePublicTool
diff --git a/InnerDetector/InDetExample/InDetRecExample/share/InDetRec_jobOptions.py b/InnerDetector/InDetExample/InDetRecExample/share/InDetRec_jobOptions.py
index 7987f21044b4cdb9b945dc5e4612f937d221445b..38558de3029b17ce906dfb2b75175ed733496a56 100755
--- a/InnerDetector/InDetExample/InDetRecExample/share/InDetRec_jobOptions.py
+++ b/InnerDetector/InDetExample/InDetRecExample/share/InDetRec_jobOptions.py
@@ -928,19 +928,20 @@ else:
# InDetTruthTrackCreation.OutputLevel = VERBOSE
topSequence += InDetTruthTrackCreation
- # --- add the truth to the truth tracks ;-)
- include ("InDetRecExample/ConfiguredInDetTrackTruth.py")
- InDetTracksTruth = ConfiguredInDetTrackTruth(InDetKeys.PseudoTracks(),
+ if InDetFlags.doSplitReco() :
+ # --- add the truth to the truth tracks ;-)
+ include ("InDetRecExample/ConfiguredInDetTrackTruth.py")
+ InDetTracksTruth = ConfiguredInDetTrackTruth(InDetKeys.PseudoTracks(),
InDetKeys.PseudoDetailedTracksTruth(),
InDetKeys.PseudoTracksTruth(),
PixelClusterTruth,
SCT_ClusterTruth,
TRT_DriftCircleTruth)
- from TrkTruthToTrack.TrkTruthToTrackConf import Trk__TruthToTrack
- InDetTruthToTrack = Trk__TruthToTrack(name = "InDetTruthToTrack",
+ from TrkTruthToTrack.TrkTruthToTrackConf import Trk__TruthToTrack
+ InDetTruthToTrack = Trk__TruthToTrack(name = "InDetTruthToTrack",
Extrapolator = TrackingCommon.getInDetExtrapolator())
- ToolSvc += InDetTruthToTrack
+ ToolSvc += InDetTruthToTrack
# Register the track collections for further processing - only if new tracking has not been running
if not InDetFlags.doNewTracking():
diff --git a/InnerDetector/InDetTruth/InDetTruthTools/src/PRD_TruthTrajectorySelectorID.cxx b/InnerDetector/InDetTruth/InDetTruthTools/src/PRD_TruthTrajectorySelectorID.cxx
index 80e7bb29a1c0a517aae60229a2d15455b44c5244..b3ab20ab83ed414a798e9c53a1a7451e46491208 100644
--- a/InnerDetector/InDetTruth/InDetTruthTools/src/PRD_TruthTrajectorySelectorID.cxx
+++ b/InnerDetector/InDetTruth/InDetTruthTools/src/PRD_TruthTrajectorySelectorID.cxx
@@ -140,16 +140,24 @@ bool InDet::PRD_TruthTrajectorySelectorID::pass( const Trk::PRD_TruthTrajectory
for ( int i = 0; i < 3 && prdIter != prdIterE; ++ prdIter ){
if( m_atlasId->is_pixel((*prdIter)->identify()) ){
const InDet::PixelCluster* pixclus=dynamic_cast<const InDet::PixelCluster*>(*prdIter);
- if (pixclus) pos.push_back (pixclus->globalPosition());
+ if (pixclus) {
+ if(!pos.empty() && (pos.back()-pixclus->globalPosition()).squaredNorm() < 9)
+ continue;
+ pos.push_back (pixclus->globalPosition());
+ }
}
else if( m_atlasId->is_sct((*prdIter)->identify()) ){
- const InDet::SCT_Cluster* sctclus=dynamic_cast<const InDet::SCT_Cluster*>(*prdIter);
- if (sctclus) pos.push_back (sctclus->globalPosition());
+ const InDet::SCT_Cluster* sctclus=dynamic_cast<const InDet::SCT_Cluster*>(*prdIter);
+ if (sctclus) {
+ if(!pos.empty() && (pos.back()-sctclus->globalPosition()).squaredNorm() < 9)
+ continue;
+ pos.push_back (sctclus->globalPosition());
+ }
}
else if( m_atlasId->is_trt((*prdIter)->identify()) ){
- continue;
+ continue;
}
- i++;
+ i++;
}
// only take trajectory if enough hits
diff --git a/Reconstruction/egamma/egammaRec/python/EMCommonRefitter.py b/Reconstruction/egamma/egammaRec/python/EMCommonRefitter.py
index da9368ef829acb9f176dcbe44ebb57fa003aeaa6..e4188ae5d8794c62b881eff6846aca05607bafad 100644
--- a/Reconstruction/egamma/egammaRec/python/EMCommonRefitter.py
+++ b/Reconstruction/egamma/egammaRec/python/EMCommonRefitter.py
@@ -32,9 +32,9 @@ def getGSFTrackFitter():
# Set up the GSF
from TrkGaussianSumFilter.TrkGaussianSumFilterConf import (
- Trk__GsfMaterialMixtureConvolution)
+ Trk__GsfMaterialMixtureConvolutionLM)
- GsfMaterialUpdator = Trk__GsfMaterialMixtureConvolution(name='GsfMaterialUpdator',MaximumNumberOfComponents=12)
+ GsfMaterialUpdator = Trk__GsfMaterialMixtureConvolutionLM(name='GsfMaterialUpdator',MaximumNumberOfComponents=12)
from TrkGaussianSumFilter.TrkGaussianSumFilterConf import (
Trk__GsfExtrapolator)
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/AlignedDynArray.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/AlignedDynArray.h
index 5809bbd07bcaac4f9a2d9c8e223f21bfb71949e8..a833fe39d9a1ccb6765c5bb640fc67ea7682567a 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/AlignedDynArray.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/AlignedDynArray.h
@@ -61,6 +61,12 @@ public:
/// size of allocated buffer
std::size_t size() const noexcept;
+
+ typedef T* iterator;
+ typedef const T* const_iterator;
+ iterator begin() { return &m_buffer[0]; }
+ iterator end() { return &m_buffer[m_size]; }
+
private:
void cleanup() noexcept;
T* m_buffer = nullptr;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/ForwardGsfFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/ForwardGsfFitter.h
index 38c10489cb58010da2a3a679283feff3a239cce1..c7c912a05008a0e92f6e259e95ede834a2e75cad 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/ForwardGsfFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/ForwardGsfFitter.h
@@ -15,6 +15,7 @@
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkFitterUtils/FitterTypes.h"
#include "TrkGaussianSumFilter/IForwardGsfFitter.h"
+#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
#include "TrkParameters/TrackParameters.h"
@@ -24,7 +25,6 @@
namespace Trk {
class IMultiStateMeasurementUpdator;
-class IMultiStateExtrapolator;
class IRIO_OnTrackCreator;
class Surface;
@@ -58,6 +58,7 @@ public:
/** Forward GSF fit using PrepRawData */
virtual std::unique_ptr<ForwardTrajectory> fitPRD(
const EventContext& ctx,
+ IMultiStateExtrapolator::Cache&,
const PrepRawDataSet&,
const TrackParameters&,
const ParticleHypothesis particleHypothesis =
@@ -66,6 +67,7 @@ public:
/** Forward GSF fit using MeasurementSet */
virtual std::unique_ptr<ForwardTrajectory> fitMeasurements(
const EventContext& ctx,
+ IMultiStateExtrapolator::Cache&,
const MeasurementSet&,
const TrackParameters&,
const ParticleHypothesis particleHypothesis =
@@ -78,6 +80,7 @@ private:
/** Progress one step along the fit */
bool stepForwardFit(
const EventContext& ctx,
+ IMultiStateExtrapolator::Cache&,
ForwardTrajectory*,
const PrepRawData*,
const MeasurementBase*,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
index 160dcea5aa5e442f31e82db531a4a802871c9989..6fc9d21ce6e9c1c2f8b08f2f2d9a419529028a4a 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
@@ -12,6 +12,7 @@
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkEventUtils/TrkParametersComparisonFunction.h"
#include "TrkFitterInterfaces/ITrackFitter.h"
+#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "TrkFitterUtils/FitterTypes.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkFitterUtils/TrackFitInputPreparator.h"
@@ -28,7 +29,6 @@
namespace Trk {
class IMultiStateMeasurementUpdator;
class MultiComponentStateOnSurface;
-class IMultiStateExtrapolator;
class IForwardGsfFitter;
class IGsfSmoother;
class FitQuality;
@@ -115,6 +115,7 @@ private:
/** Produces a perigee from a smoothed trajectory */
const MultiComponentStateOnSurface* makePerigee(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache&,
const SmoothedTrajectory*,
const ParticleHypothesis particleHypothesis = nonInteracting) const;
@@ -150,7 +151,14 @@ private:
"GsfSmoother",
"Trk::GsfSmoother/GsfSmoother",
"" };
-
+
+ Gaudi::Property<bool> m_StoreMCSOS{
+ this,
+ "StoreMCSOS",
+ false,
+ "Store multicomponent state or single state in final trajectory"
+ };
+
bool m_reintegrateOutliers;
bool m_makePerigee;
bool m_refitOnMeasurementBase;
@@ -173,6 +181,10 @@ private:
mutable std::atomic<int> m_PerigeeFailure;
// Number of Tracks that fail fit Quailty test
mutable std::atomic<int> m_fitQualityFailure;
+ // Number of Tracks that are successfull
+ mutable std::atomic<int> m_fitSuccess;
+
+
};
} // end Trk namespace
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfExtrapolator.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfExtrapolator.h
index 350e84992f3e16b8bc0299d7cf7a4b39ea44fc84..8db31e2430a36034a5fa861ef841bc2f2ec40fa9 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfExtrapolator.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfExtrapolator.h
@@ -19,6 +19,7 @@
#include "TrkExInterfaces/INavigator.h"
#include "TrkExInterfaces/IPropagator.h"
#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
+#include "TrkGaussianSumFilter/IMultiStateMaterialEffects.h"
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/IChronoStatSvc.h"
@@ -45,37 +46,6 @@ class MaterialProperties;
class IMultiComponentStateMerger;
class IMaterialMixtureConvolution;
class IMultipleScatteringUpdator;
-/** @struct StateAtBoundarySurface
- - Structure to contain information about a state at the interface between
- tracking volumes
- */
-
-struct StateAtBoundarySurface
-{
-
- /** Data members */
- const MultiComponentState* stateAtBoundary;
- const TrackParameters* navigationParameters;
- const TrackingVolume* trackingVolume;
-
- /** Default constructor */
- StateAtBoundarySurface()
- : stateAtBoundary(nullptr)
- , navigationParameters(nullptr)
- , trackingVolume(nullptr)
- {}
-
- /** Update State at Boundary Surface Information */
- void updateBoundaryInformation(const MultiComponentState* boundaryState,
- const TrackParameters* navParameters,
- const TrackingVolume* nextVolume)
- {
- stateAtBoundary = boundaryState;
- navigationParameters = navParameters;
- trackingVolume = nextVolume;
- }
-};
-
/** @class GsfExtrapolator */
class GsfExtrapolator
@@ -99,6 +69,7 @@ public:
/** Configured AlgTool extrapolation method (1) */
virtual MultiComponentState extrapolate(
const EventContext& ctx,
+ Cache&,
const MultiComponentState&,
const Surface&,
PropDirection direction = anyDirection,
@@ -125,33 +96,6 @@ public:
ParticleHypothesis particle) const override final;
private:
- struct Cache
- {
- bool m_recall; //!< Flag the recall solution
- const Surface* m_recallSurface; //!< Surface for recall
- const Layer* m_recallLayer; //!< Layer for recall
- const TrackingVolume*
- m_recallTrackingVolume; //!< Tracking volume for recall
- StateAtBoundarySurface
- m_stateAtBoundarySurface; //!< Instance of structure describing the state
- //!< at a boundary of tracking volumes
- std::unique_ptr<std::vector<const Trk::TrackStateOnSurface*>> m_matstates;
- std::vector<std::unique_ptr<const MultiComponentState>>
- m_mcsGarbageBin; //!< Garbage bin for MultiComponentState objects
- std::vector<std::unique_ptr<const TrackParameters>>
- m_tpGarbageBin; //!< Garbage bin for TrackParameter objects
-
- Cache()
- : m_recall(false)
- , m_recallSurface(nullptr)
- , m_recallLayer(nullptr)
- , m_recallTrackingVolume(nullptr)
- , m_stateAtBoundarySurface()
- , m_matstates(nullptr)
- , m_mcsGarbageBin()
- , m_tpGarbageBin()
- {}
- };
/** These are the methods that do the actual heavy lifting when extrapolating
* with a cache */
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolution.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolution.h
index 0c32a8a0c66221808eba97fb06c2371d76ef2415..65d4e5296cd14e1e7615df4077ffd23404a5b792 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolution.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolution.h
@@ -44,21 +44,24 @@ public:
virtual StatusCode finalize() override;
//!< Convolution with full material properties
- virtual MultiComponentState update(const MultiComponentState&,
+ virtual MultiComponentState update(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
ParticleHypothesis particleHypothesis =
nonInteracting) const override final;
//!< Convolution with pre-measurement-update material properties
- virtual MultiComponentState preUpdate(const MultiComponentState&,
+ virtual MultiComponentState preUpdate(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
ParticleHypothesis particleHypothesis =
nonInteracting) const override final;
//!< Convolution with post-measurement-update material properties
- virtual MultiComponentState postUpdate(const MultiComponentState&,
+ virtual MultiComponentState postUpdate(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
ParticleHypothesis particleHypothesis =
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolutionLM.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolutionLM.h
new file mode 100644
index 0000000000000000000000000000000000000000..4611856669cdf962b39a5f1f664ba9b1fa9037d5
--- /dev/null
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfMaterialMixtureConvolutionLM.h
@@ -0,0 +1,138 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * @file GsfMaterialMixtureConvolutionLM.h
+ * @date Thursday 7th September 2006
+ * @author Tom Athkinson, Anthony Morley, Christos Anastopoulos
+ * @brief Class description for convolution of GSF material mixture
+ */
+
+#ifndef TrkGsfMaterialMixtureConvolutionLM_H
+#define TrkGsfMaterialMixtureConvolutionLM_H
+
+#include "TrkGaussianSumFilter/IMaterialMixtureConvolution.h"
+#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+namespace Trk {
+
+class IMultiStateMaterialEffectsUpdator;
+class IMultiStateMaterialEffects;
+class Layer;
+class MaterialProperties;
+
+class GsfMaterialMixtureConvolutionLM
+ : public AthAlgTool
+ , virtual public IMaterialMixtureConvolution
+{
+
+public:
+
+ enum MaterialUpdateType{ Normal = 0, Preupdate=1, Postupdate=2 };
+
+ //!< Constructor with AlgTool parameters
+ GsfMaterialMixtureConvolutionLM(const std::string&,
+ const std::string&,
+ const IInterface*);
+
+ //!< Destructor
+ virtual ~GsfMaterialMixtureConvolutionLM();
+
+ //!< AlgTool initialise method
+ virtual StatusCode initialize() override;
+
+ //!< AlgTool finalize method
+ virtual StatusCode finalize() override;
+
+ //!< Convolution with full material properties
+ virtual MultiComponentState update(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
+ const Layer&,
+ PropDirection direction = anyDirection,
+ ParticleHypothesis particleHypothesis =
+ nonInteracting) const override final;
+
+ //!< Convolution with pre-measurement-update material properties
+ virtual MultiComponentState preUpdate(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
+ const Layer&,
+ PropDirection direction = anyDirection,
+ ParticleHypothesis particleHypothesis =
+ nonInteracting) const override final;
+
+ //!< Convolution with post-measurement-update material properties
+ virtual MultiComponentState postUpdate(std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const MultiComponentState&,
+ const Layer&,
+ PropDirection direction = anyDirection,
+ ParticleHypothesis particleHypothesis =
+ nonInteracting) const override final;
+
+ //!< Retain for now redundant simplified material effects
+ virtual MultiComponentState simplifiedMaterialUpdate(
+ const MultiComponentState& ,
+ PropDirection,
+ ParticleHypothesis) const override final { return {}; };
+
+
+ private:
+ Gaudi::Property<unsigned int> m_maximumNumberOfComponents{
+ this,
+ "MaximumNumberOfComponents",
+ 12,
+ "Maximum number of components"
+ };
+
+ ToolHandle<IMultiStateMaterialEffectsUpdator> m_updator{
+ this,
+ "MaterialEffectsUpdator",
+ "Trk::GsfMaterialEffectsUpdator/GsfMaterialEffectsUpdator",
+ ""
+ };
+
+ ToolHandle<IMultiStateMaterialEffects> m_materialEffects{
+ this,
+ "MaterialEffects",
+ "Trk::GsfCombinedMaterialEffects/GsfCombinedMaterialEffects",
+ ""
+ };
+
+ Gaudi::Property<bool> m_useReferenceMaterial{
+ this,
+ "UseReferenceMaterial",
+ false,
+ ""
+ };
+
+ Gaudi::Property<double> m_momentumCut{
+ this,
+ "MinimalMomentum",
+ 250. * Gaudi::Units::MeV,
+ ""
+ };
+
+
+ Trk::MultiComponentState update(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
+ const Trk::MultiComponentState& inputState,
+ const Trk::Layer& layer,
+ Trk::PropDirection direction,
+ Trk::ParticleHypothesis particleHypothesis,
+ MaterialUpdateType updateType) const;
+
+ bool updateP(double& qOverP, double deltaP) const;
+
+ std::pair< const Trk::MaterialProperties*, double > getMaterialProperties(
+ const Trk::TrackParameters* trackParameters,
+ const Trk::Layer& layer ) const;
+
+
+};
+
+} // end Trk namespace
+
+#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfSmoother.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfSmoother.h
index f7a363901c464eb4d9df90c6bbe09791d36783ac..615054c82b7fbf2d162d313f26cd53d1bc6c6e7c 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfSmoother.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfSmoother.h
@@ -54,6 +54,7 @@ public:
/** Gsf smoother method */
virtual SmoothedTrajectory* fit(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache&,
const ForwardTrajectory&,
const ParticleHypothesis particleHypothesis = nonInteracting,
const CaloCluster_OnTrack* ccot = nullptr) const;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IForwardGsfFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IForwardGsfFitter.h
index ebc4bbd4c79144dd3670c40c46a966ac6ab82965..f514deb5402bd767c625b094b299494007814235 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IForwardGsfFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IForwardGsfFitter.h
@@ -18,6 +18,7 @@
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkFitterUtils/FitterTypes.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "TrkParameters/TrackParameters.h"
#include "GaudiKernel/EventContext.h"
@@ -25,7 +26,6 @@
namespace Trk {
class IMultiStateMeasurementUpdator;
-class IMultiStateExtrapolator;
class IRIO_OnTrackCreator;
class Surface;
@@ -56,6 +56,7 @@ public:
/** Forward GSF fit using PrepRawData */
virtual std::unique_ptr<ForwardTrajectory> fitPRD(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache&,
const PrepRawDataSet&,
const TrackParameters&,
const ParticleHypothesis particleHypothesis = nonInteracting) const = 0;
@@ -63,6 +64,7 @@ public:
/** Forward GSF fit using MeasurementSet */
virtual std::unique_ptr<ForwardTrajectory> fitMeasurements(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache&,
const MeasurementSet&,
const TrackParameters&,
const ParticleHypothesis particleHypothesis = nonInteracting) const = 0;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfSmoother.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfSmoother.h
index 2dbf97f0af0fd275c95e24983298a5174fb6d4ec..11ae1e16db27e39c95f90bd4cf7c9a80bd2643c0 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfSmoother.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfSmoother.h
@@ -15,6 +15,7 @@
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkFitterUtils/FitterTypes.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "GaudiKernel/EventContext.h"
#include "GaudiKernel/IAlgTool.h"
@@ -23,7 +24,6 @@
namespace Trk {
class IMultiStateMeasurementUpdator;
-class IMultiStateExtrapolator;
class CaloCluster_OnTrack;
static const InterfaceID InterfaceID_GsfSmoother("GsfSmoother", 1, 0);
@@ -48,6 +48,7 @@ public:
/** Gsf smoother method */
virtual SmoothedTrajectory* fit(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache&,
const ForwardTrajectory&,
const ParticleHypothesis particleHypothesis = nonInteracting,
const CaloCluster_OnTrack* ccot = nullptr) const = 0;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMaterialMixtureConvolution.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMaterialMixtureConvolution.h
index 10aa90246091ea06c7bc9e960995c57c660732c0..469b19c66019630148af1f0d72a00cad1739350b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMaterialMixtureConvolution.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMaterialMixtureConvolution.h
@@ -16,6 +16,7 @@
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+#include "TrkGaussianSumFilter/IMultiStateMaterialEffects.h"
namespace Trk {
class Layer;
@@ -37,6 +38,7 @@ public:
//!< Convolution with full material properties
virtual MultiComponentState update(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
@@ -44,6 +46,7 @@ public:
//!< Convolution with pre-measurement-update material properties
virtual MultiComponentState preUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
@@ -51,6 +54,7 @@ public:
//!< Convolution with post-measurement-update material properties
virtual MultiComponentState postUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const MultiComponentState&,
const Layer&,
PropDirection direction = anyDirection,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateExtrapolator.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateExtrapolator.h
index 42d6af11468a19e08fe38af862c45f7d12be79b3..d7e0a75b14615ed3233ec68643e91279fddb49f9 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateExtrapolator.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateExtrapolator.h
@@ -14,6 +14,7 @@
#define TrkIMultiStateExtrapolator_H
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+#include "TrkGaussianSumFilter/IMultiStateMaterialEffects.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkEventPrimitives/PropDirection.h"
@@ -35,6 +36,32 @@ class Track;
class TrackingVolume;
class TrackStateOnSurface;
+/** @struct StateAtBoundarySurface
+ - Structure to contain information about a state at the interface between
+ tracking volumes
+ */
+
+struct StateAtBoundarySurface
+{
+
+ /** Data members */
+ const MultiComponentState* stateAtBoundary = nullptr;
+ const TrackParameters* navigationParameters = nullptr;
+ const TrackingVolume* trackingVolume = nullptr;
+
+
+ /** Update State at Boundary Surface Information */
+ void updateBoundaryInformation(const MultiComponentState* boundaryState,
+ const TrackParameters* navParameters,
+ const TrackingVolume* nextVolume)
+ {
+ stateAtBoundary = boundaryState;
+ navigationParameters = navParameters;
+ trackingVolume = nextVolume;
+ }
+};
+
+
static const InterfaceID IID_IMultiStateExtrapolator("IMultiStateExtrapolator",
1,
0);
@@ -42,7 +69,52 @@ static const InterfaceID IID_IMultiStateExtrapolator("IMultiStateExtrapolator",
class IMultiStateExtrapolator : virtual public IAlgTool
{
public:
- /** Virtual destructor */
+ /** @brief MultiStateExtrapolator cache class.
+ * This object holds information regarding the state of the extrpolation process
+ * as well as a large store for material effects properties.
+ *
+ * The is owned by the tools that call the extrapolatator (Currently the GSFSmoother or Forward fitter)
+ * and it not exposed to the caller of the track fitter.
+ * */
+ struct Cache
+ {
+ bool m_recall = false; //!< Flag the recall solution
+ const Surface* m_recallSurface = nullptr; //!< Surface for recall
+ const Layer* m_recallLayer = nullptr; //!< Layer for recall
+ const TrackingVolume*
+ m_recallTrackingVolume = nullptr; //!< Tracking volume for recall
+ StateAtBoundarySurface
+ m_stateAtBoundarySurface; //!< Instance of structure describing the state
+ //!< at a boundary of tracking volumes
+ std::unique_ptr<std::vector<const Trk::TrackStateOnSurface*>> m_matstates;
+ std::vector<std::unique_ptr<const MultiComponentState>>
+ m_mcsGarbageBin; //!< Garbage bin for MultiComponentState objects
+ std::vector<std::unique_ptr<const TrackParameters>>
+ m_tpGarbageBin; //!< Garbage bin for TrackParameter objects
+
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>
+ m_materialEffectsCaches; //!< Large cache for material effects
+
+
+ Cache()
+ {
+ m_materialEffectsCaches.reserve(72);
+ }
+
+ void reset(){
+ m_recall = false;
+ m_recallSurface = nullptr;
+ m_recallLayer = nullptr;
+ m_recallTrackingVolume = nullptr;
+ m_matstates.reset(nullptr);
+ m_stateAtBoundarySurface.updateBoundaryInformation(nullptr,nullptr,nullptr);
+ m_mcsGarbageBin.clear();
+ m_tpGarbageBin.clear();
+ m_materialEffectsCaches.clear();
+ };
+ };
+
+ /** Virtual destructor */
virtual ~IMultiStateExtrapolator() = default;
/** AlgTool interface method */
@@ -54,6 +126,7 @@ public:
/** Configured AlgTool extrapolation method (1) */
virtual MultiComponentState extrapolate(
const EventContext& ctx,
+ Cache&,
const MultiComponentState&,
const Surface&,
PropDirection direction = anyDirection,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateMaterialEffects.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateMaterialEffects.h
index 1d235b978b8dc370681aace97f95aad0e539e2e6..870ee67246dd40e7cba3c32a92608f3de439fb6b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateMaterialEffects.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiStateMaterialEffects.h
@@ -44,6 +44,7 @@ public:
{
weights.reserve(6);
deltaPs.reserve(6);
+ deltaParameters.reserve(6);
deltaCovariances.reserve(6);
}
Cache(Cache&&) = default;
@@ -60,14 +61,29 @@ public:
* "you must use the Eigen::aligned_allocator (not another aligned
* allocator), and #include <Eigen/StdVector>."
*/
+ std::vector<AmgVector(5),Eigen::aligned_allocator<AmgVector(5)>>
+ deltaParameters;
std::vector<AmgSymMatrix(5), Eigen::aligned_allocator<AmgSymMatrix(5)>>
deltaCovariances;
void reset()
{
weights.clear();
deltaPs.clear();
+ deltaParameters.clear();
deltaCovariances.clear();
}
+ void resetAndAddDummyValues()
+ {
+ reset();
+ weights.push_back(1);
+ deltaPs.push_back(0);
+ AmgVector(5) newParameters;
+ newParameters.setZero();
+ deltaParameters.push_back( std::move(newParameters) );
+ AmgSymMatrix(5) newCovarianceMatrix;
+ newCovarianceMatrix.setZero();
+ deltaCovariances.push_back( std::move(newCovarianceMatrix) );
+ }
};
/** Alg tool and IAlgTool interface method */
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
index d18a5aee90fa99f9e93538c75bb9cb0f31a6a7ad..95cd0d23d8a861652448da81a774ab4350f29a9d 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
@@ -31,6 +31,26 @@ void
combineWithWeight(Trk::ComponentParameters& mergeTo,
const Trk::ComponentParameters& addThis);
+/** @brief Combined/merge a component to the parts of another */
+void
+combineWithWeight( Trk::ComponentParameters& mergeTo,
+ const AmgVector(5)& secondParameters,
+ const AmgSymMatrix(5)* secondMeasuredCov,
+ const double secondWeight );
+
+
+/** @brief Update first parameters */
+void combineParametersWithWeight(
+ AmgVector(5)& firstParameters, double& firstWeight,
+ const AmgVector(5)& secondParameters, const double secondWeight );
+
+/** @brief Update cov matrix */
+void combineCovWithWeight(
+ const AmgVector(5)& firstParameters, AmgSymMatrix(5)& firstMeasuredCov, const double firstWeight,
+ const AmgVector(5)& secondParameters, const AmgSymMatrix(5)& secondMeasuredCov, const double secondWeight );
+
+
+
/** @brief Calculate combined state and weight of many components */
std::unique_ptr<Trk::ComponentParameters>
combineWithWeight(const MultiComponentState&,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/ForwardGsfFitter.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/ForwardGsfFitter.cxx
index 12ab02457be411564eac4fe1300c78af4ee5336b..f317764a29d4be26782b557a60498cc62d6948e0 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/ForwardGsfFitter.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/ForwardGsfFitter.cxx
@@ -91,6 +91,7 @@ Trk::ForwardGsfFitter::configureTools(
std::unique_ptr<Trk::ForwardTrajectory>
Trk::ForwardGsfFitter::fitPRD(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
const Trk::PrepRawDataSet& inputPrepRawDataSet,
const Trk::TrackParameters& estimatedTrackParametersNearOrigin,
const Trk::ParticleHypothesis particleHypothesis) const
@@ -149,6 +150,7 @@ Trk::ForwardGsfFitter::fitPRD(
// stepForwardFit method is updated
bool stepIsValid = stepForwardFit(
ctx,
+ extrapolatorCache,
forwardTrajectory.get(),
*prepRawData,
nullptr,
@@ -172,6 +174,7 @@ Trk::ForwardGsfFitter::fitPRD(
std::unique_ptr<Trk::ForwardTrajectory>
Trk::ForwardGsfFitter::fitMeasurements(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
const Trk::MeasurementSet& inputMeasurementSet,
const Trk::TrackParameters& estimatedTrackParametersNearOrigin,
const Trk::ParticleHypothesis particleHypothesis) const
@@ -230,6 +233,7 @@ Trk::ForwardGsfFitter::fitMeasurements(
for (; measurement != inputMeasurementSet.end(); ++measurement) {
bool stepIsValid = stepForwardFit(ctx,
+ extrapolatorCache,
forwardTrajectory.get(),
nullptr,
*measurement,
@@ -252,6 +256,7 @@ Trk::ForwardGsfFitter::fitMeasurements(
bool
Trk::ForwardGsfFitter::stepForwardFit(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
ForwardTrajectory* forwardTrajectory,
const Trk::PrepRawData* originalPrepRawData,
const Trk::MeasurementBase* originalMeasurement,
@@ -276,6 +281,7 @@ Trk::ForwardGsfFitter::stepForwardFit(
// =================================================================
Trk::MultiComponentState extrapolatedState =
m_extrapolator->extrapolate(ctx,
+ extrapolatorCache,
updatedState,
surface,
Trk::alongMomentum,
@@ -284,6 +290,8 @@ Trk::ForwardGsfFitter::stepForwardFit(
if (extrapolatedState.empty()) {
ATH_MSG_DEBUG("Extrapolation failed... returning false");
return false;
+ } else {
+ ATH_MSG_DEBUG("Extrapolation worked... state size: "<< extrapolatedState.size() );
}
// =======================
// Measurement Preparation
@@ -321,6 +329,8 @@ Trk::ForwardGsfFitter::stepForwardFit(
if (updatedState.empty()) {
ATH_MSG_DEBUG("Measurement update of the state failed... Exiting!");
return false;
+ } else {
+ ATH_MSG_DEBUG("Measurement update of the state worked : " << updatedState.size() );
}
// Bail if the fit quality is not defined:
if (!fitQuality) {
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
index 618c26e733a173422b4b9232762e0a1275b14aa1..b2654163bd4876b21e864738fd70861272ddb338 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
@@ -55,6 +55,7 @@ Trk::GaussianSumFitter::GaussianSumFitter(const std::string& type,
, m_SmootherFailure{ 0 }
, m_PerigeeFailure{ 0 }
, m_fitQualityFailure{ 0 }
+ , m_fitSuccess{ 0 }
{
declareInterface<ITrackFitter>(this);
@@ -162,6 +163,8 @@ Trk::GaussianSumFitter::finalize()
<< "Number of MakePerigee Failures: " << m_PerigeeFailure << '\n'
<< "Number of Trks that fail fitquality test: " << m_fitQualityFailure
<< '\n'
+ << "Number of successful fits: " << m_fitSuccess << '\n'
+ << '\n'
<< "-----------------------------------------------" << '\n'
<< "Finalisation of " << name() << " was successful");
return StatusCode::SUCCESS;
@@ -323,11 +326,15 @@ Trk::GaussianSumFitter::fit(
sortedPrepRawDataSet.end(),
prdComparisonFunction);
}
+ // Create Extrapolator cache that holds material effects cache;
+ Trk::IMultiStateExtrapolator::Cache extrapolatorCache;
+
// Perform GSF forwards fit
ForwardTrajectory* forwardTrajectory =
m_forwardGsfFitter
->fitPRD(ctx,
+ extrapolatorCache,
sortedPrepRawDataSet,
estimatedParametersNearOrigin,
particleHypothesis)
@@ -348,7 +355,7 @@ Trk::GaussianSumFitter::fit(
// Perform GSF smoother operation
SmoothedTrajectory* smoothedTrajectory =
- m_gsfSmoother->fit(ctx,*forwardTrajectory, particleHypothesis);
+ m_gsfSmoother->fit(ctx, extrapolatorCache, *forwardTrajectory, particleHypothesis);
// Protect against failed smoother fit
if (!smoothedTrajectory) {
@@ -370,7 +377,7 @@ Trk::GaussianSumFitter::fit(
if (m_makePerigee) {
const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface =
- this->makePerigee(ctx,smoothedTrajectory, particleHypothesis);
+ makePerigee(ctx, extrapolatorCache, smoothedTrajectory, particleHypothesis);
ATH_MSG_DEBUG(
"perigeeMultiStateOnSurface :" << perigeeMultiStateOnSurface);
if (perigeeMultiStateOnSurface) {
@@ -389,11 +396,28 @@ Trk::GaussianSumFitter::fit(
delete forwardTrajectory;
// Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
+
+ // Store only TSOS in tracks instead of MCSOS
+ if( !m_StoreMCSOS ){
+ auto slimmedSmoothedTrajectory = std::make_unique<Trk::SmoothedTrajectory>();
+ for( const Trk::TrackStateOnSurface* tsos : *smoothedTrajectory ){
+ slimmedSmoothedTrajectory->push_back( new Trk::TrackStateOnSurface( *tsos ) );
+ }
+ delete smoothedTrajectory;
+ // Create new track
+ Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
+ info.setTrackProperties(TrackInfo::BremFit);
+ info.setTrackProperties(TrackInfo::BremFitSuccessful);
+ ++m_fitSuccess;
+ return std::make_unique<Track>(info, slimmedSmoothedTrajectory.release(), fitQuality);
+ }
+
// Create new track
Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
info.setTrackProperties(TrackInfo::BremFit);
info.setTrackProperties(TrackInfo::BremFitSuccessful);
+ ++m_fitSuccess;
return std::make_unique<Track>(info, smoothedTrajectory, fitQuality);
}
@@ -465,10 +489,17 @@ Trk::GaussianSumFitter::fit(
measurementBaseComparisonFunction);
}
+
+ // Create Extrapolator cache that holds material effects cache;
+ Trk::IMultiStateExtrapolator::Cache extrapolatorCache;
+
+
+
// Perform GSF forwards fit - new memory allocated in forwards fitter
ForwardTrajectory* forwardTrajectory =
m_forwardGsfFitter
->fitMeasurements(ctx,
+ extrapolatorCache,
sortedMeasurementSet,
estimatedParametersNearOrigin,
particleHypothesis)
@@ -490,7 +521,7 @@ Trk::GaussianSumFitter::fit(
// Perform GSF smoother operation
SmoothedTrajectory* smoothedTrajectory =
- m_gsfSmoother->fit(ctx, *forwardTrajectory, particleHypothesis, ccot);
+ m_gsfSmoother->fit(ctx, extrapolatorCache, *forwardTrajectory, particleHypothesis, ccot);
// Protect against failed smoother fit
if (!smoothedTrajectory) {
@@ -513,7 +544,7 @@ Trk::GaussianSumFitter::fit(
if (m_makePerigee) {
const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface =
- this->makePerigee(ctx,smoothedTrajectory, particleHypothesis);
+ makePerigee(ctx, extrapolatorCache, smoothedTrajectory, particleHypothesis);
ATH_MSG_DEBUG(
"perigeeMultiStateOnSurface :" << perigeeMultiStateOnSurface);
@@ -535,10 +566,27 @@ Trk::GaussianSumFitter::fit(
// Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
+
+ // Store only TSOS in tracks instead of MCSOS
+ if( !m_StoreMCSOS ){
+ auto slimmedSmoothedTrajectory = std::make_unique<Trk::SmoothedTrajectory>();
+ for( const Trk::TrackStateOnSurface* tsos : *smoothedTrajectory ){
+ slimmedSmoothedTrajectory->push_back( new Trk::TrackStateOnSurface(*tsos) );
+ }
+ delete smoothedTrajectory;
+ // Create new track
+ Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
+ info.setTrackProperties(TrackInfo::BremFit);
+ info.setTrackProperties(TrackInfo::BremFitSuccessful);
+ ++m_fitSuccess;
+ return std::make_unique<Track>(info, slimmedSmoothedTrajectory.release(), fitQuality);
+ }
+
// Create new track
Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
info.setTrackProperties(TrackInfo::BremFit);
info.setTrackProperties(TrackInfo::BremFitSuccessful);
+ ++m_fitSuccess;
return std::make_unique<Track>(info, smoothedTrajectory, fitQuality);
}
@@ -698,6 +746,7 @@ Trk::GaussianSumFitter::fit(const EventContext& ctx,
const Trk::MultiComponentStateOnSurface*
Trk::GaussianSumFitter::makePerigee(
const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
const Trk::SmoothedTrajectory* smoothedTrajectory,
const Trk::ParticleHypothesis particleHypothesis) const
{
@@ -733,6 +782,7 @@ Trk::GaussianSumFitter::makePerigee(
// Extrapolate to perigee, taking material effects considerations into account
Trk::MultiComponentState stateExtrapolatedToPerigee =
m_extrapolator->extrapolate(ctx,
+ extrapolatorCache,
*multiComponentState,
perigeeSurface,
m_directionToPerigee,
@@ -763,7 +813,7 @@ Trk::GaussianSumFitter::makePerigee(
pattern(0);
pattern.set(Trk::TrackStateOnSurface::Perigee);
- if (fabs(combinedPerigee->parameters()[Trk::qOverP]) > 1e8) {
+ if (std::abs(combinedPerigee->parameters()[Trk::qOverP]) > 1e8) {
// Protection against 0-momentum track .. this check should NEVER be needed.
//actual cutoff is 0.01eV track
ATH_MSG_ERROR(
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
index c5b2735ab7536d41f3b4db68959c1bbb94b2ade6..b53ae153a282b0f8ff72761696eaaa9746eb3f2d 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
@@ -72,7 +72,7 @@ Trk::GsfCombinedMaterialEffects::compute(
Trk::ParticleHypothesis particleHypothesis) const
{
- ATH_MSG_DEBUG("Computing combined material effects");
+ ATH_MSG_DEBUG("Computing combined material effects, P : " << componentParameters.first->momentum().norm() << " W " << componentParameters.second );
// Reset everything before computation
cache.reset();
@@ -163,7 +163,6 @@ Trk::GsfCombinedMaterialEffects::compute(
(*multipleScatter_weightsIterator) * (*energyLoss_weightsIterator);
double combinedDeltaP =
(*multipleScatter_deltaPsIterator) + (*energyLoss_deltaPsIterator);
-
cache.weights.push_back(combinedWeight);
cache.deltaPs.push_back(combinedDeltaP);
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfExtrapolator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfExtrapolator.cxx
index 0dfb08aefafa805c97f66a89de89ec7343c3b7e0..251154f86a8996fae6eb3031eb7ece63562786c1 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfExtrapolator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfExtrapolator.cxx
@@ -145,6 +145,7 @@ Trk::GsfExtrapolator::extrapolateImpl(
const Trk::BoundaryCheck& boundaryCheck,
Trk::ParticleHypothesis particleHypothesis) const
{
+ ATH_MSG_DEBUG("Calling extrpolate: " << multiComponentState.size() );
auto buff_extrapolateCalls = m_extrapolateCalls.buffer();
// If the extrapolation is to be without material effects simply revert to the
@@ -312,7 +313,7 @@ Trk::GsfExtrapolator::extrapolateImpl(
double revisedDistance =
(referenceParameters->position() - newDestination).mag();
- double distanceChange = fabs(revisedDistance - initialDistance);
+ double distanceChange = std::abs(revisedDistance - initialDistance);
if (revisedDistance > initialDistance && distanceChange > 0.01) {
ATH_MSG_DEBUG("Navigation break. Initial separation: "
@@ -432,6 +433,7 @@ Trk::GsfExtrapolator::extrapolateImpl(
return {};
}
// After successful extrapolation return the state
+ ATH_MSG_DEBUG("extrapolateInsideVolume() successful: " << destinationState.size() );
return destinationState;
}
/*
@@ -469,17 +471,19 @@ Trk::GsfExtrapolator::extrapolateDirectlyImpl(
Trk::MultiComponentState
Trk::GsfExtrapolator::extrapolate(
const EventContext& ctx,
+ Cache& cache,
const Trk::MultiComponentState& multiComponentState,
const Trk::Surface& surface,
Trk::PropDirection direction,
const Trk::BoundaryCheck& boundaryCheck,
Trk::ParticleHypothesis particleHypothesis) const
{
- Cache cache{};
if (multiComponentState.empty()) {
ATH_MSG_DEBUG("MultiComponentState is empty...");
return {};
}
+
+ cache.reset();
// Set the propagator to that one corresponding to the configuration level
const Trk::IPropagator* currentPropagator =
@@ -629,7 +633,7 @@ Trk::GsfExtrapolator::extrapolateToVolumeBoundary(
else if (trackingVolume.confinedLayers() &&
associatedLayer->layerMaterialProperties()) {
Trk::MultiComponentState updatedState = m_materialUpdator->postUpdate(
- *currentState, *layer, direction, particleHypothesis);
+ cache.m_materialEffectsCaches, *currentState, *layer, direction, particleHypothesis);
if (!updatedState.empty()) {
addMaterialtoVector(cache, layer, currentState->begin()->first.get());
@@ -752,7 +756,7 @@ Trk::GsfExtrapolator::extrapolateToVolumeBoundary(
ATH_MSG_DEBUG("Boundary surface has material - updating properties");
assert(currentState);
matUpdatedState = m_materialUpdator->postUpdate(
- *currentState, *layerAtBoundary, direction, particleHypothesis);
+ cache.m_materialEffectsCaches, *currentState, *layerAtBoundary, direction, particleHypothesis);
}
}
@@ -809,7 +813,7 @@ Trk::GsfExtrapolator::extrapolateInsideVolume(
{
ATH_MSG_DEBUG("GSF extrapolateInsideVolume() in tracking volume: "
- << trackingVolume.volumeName());
+ << trackingVolume.volumeName() << " with " << multiComponentState.size() << " components" );
/*
* We use current State to track where we are
@@ -856,7 +860,7 @@ Trk::GsfExtrapolator::extrapolateInsideVolume(
associatedLayer->layerMaterialProperties()) {
updatedState = m_materialUpdator->postUpdate(
- *currentState, *associatedLayer, direction, particleHypothesis);
+ cache.m_materialEffectsCaches, *currentState, *associatedLayer, direction, particleHypothesis);
if (!updatedState.empty()) {
addMaterialtoVector(
@@ -1090,7 +1094,7 @@ Trk::GsfExtrapolator::extrapolateToIntermediateLayer(
------------------------------------- */
Trk::MultiComponentState updatedState = m_materialUpdator->update(
- destinationState, layer, direction, particleHypothesis);
+ cache.m_materialEffectsCaches, destinationState, layer, direction, particleHypothesis);
if (updatedState.empty()) {
return destinationState;
@@ -1166,9 +1170,11 @@ Trk::GsfExtrapolator::extrapolateToDestinationLayer(
Trk::MultiComponentState updatedState{};
if (startLayer != &layer) {
updatedState = m_materialUpdator->preUpdate(
- destinationState, layer, direction, particleHypothesis);
+ cache.m_materialEffectsCaches, destinationState, layer, direction, particleHypothesis);
}
+ ATH_MSG_DEBUG( "State size after preUpdate: " << updatedState.size() );
+
if (updatedState.empty()) {
return destinationState;
}
@@ -1201,7 +1207,7 @@ Trk::GsfExtrapolator::extrapolateSurfaceBasedMaterialEffects(
// Check the multi-component state is populated
if (multiComponentState.empty()) {
- ATH_MSG_WARNING("Multi component state passed to extrapolateInsideVolume "
+ ATH_MSG_WARNING("Multi component state passed to extrapolateSurfaceBasedMaterialEffects "
"is not populated... returning 0");
return {};
}
@@ -1288,7 +1294,7 @@ Trk::GsfExtrapolator::multiStatePropagate(
}
ATH_MSG_DEBUG("GSF multiStatePropagate() propagated "
- << propagatedState.size() << "components");
+ << propagatedState.size() << " components");
// Protect against empty propagation
if (propagatedState.empty() || sumw < 0.1) {
ATH_MSG_DEBUG("multiStatePropagate failed... ");
@@ -1320,7 +1326,7 @@ Trk::GsfExtrapolator::propagatorType(
// field and material properties ST : the following check may fail as the dEdX
// is often dummy for dense volumes - switch to rho or zOverAtimesRho ?
unsigned int propagatorMode =
- (magneticFieldMode > 1 && fabs(trackingVolume.dEdX) < 10e-2) ? 2 : 3;
+ (magneticFieldMode > 1 && std::abs(trackingVolume.dEdX) < 10e-2) ? 2 : 3;
unsigned int returnType = (propagatorMode > m_propagatorConfigurationLevel)
? m_propagatorConfigurationLevel
@@ -1352,7 +1358,7 @@ Trk::GsfExtrapolator::initialiseNavigation(
Trk::PropDirection direction) const
{
- ATH_MSG_DEBUG("initialiseNavigation !!!");
+ ATH_MSG_DEBUG("initialiseNavigation !!! : " << multiComponentState.size() );
// Empty the garbage bin
ATH_MSG_DEBUG("Destination to surface [r,z] ["
<< surface.center().perp() << ",\t" << surface.center().z()
@@ -1505,7 +1511,7 @@ Trk::GsfExtrapolator::addMaterialtoVector(Cache& cache,
// Determine the pathCorrection if the material properties exist
pathcorr = materialProperties
- ? 1. / fabs(surface->normal().dot(nextPar->momentum().unit()))
+ ? 1. / std::abs(surface->normal().dot(nextPar->momentum().unit()))
: 0.;
}
@@ -1527,7 +1533,7 @@ Trk::GsfExtrapolator::addMaterialtoVector(Cache& cache,
nextPar->position(), nextPar->momentum());
double thick = pathcorr * materialProperties->thickness();
double dInX0 = thick / materialProperties->x0();
- double absP = 1 / fabs(nextPar->parameters()[Trk::qOverP]);
+ double absP = 1 / std::abs(nextPar->parameters()[Trk::qOverP]);
double scatsigma = sqrt(
m_msupdators->sigmaSquare(*materialProperties, absP, pathcorr, particle));
Trk::ScatteringAngles* newsa = new Trk::ScatteringAngles(
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialEffectsUpdator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialEffectsUpdator.cxx
index 9657a682b8513f15c507950d39186a9ba9f1e9f8..16c18501e86bf97244cb50a0a593fc66e978b04a 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialEffectsUpdator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialEffectsUpdator.cxx
@@ -98,7 +98,7 @@ Trk::GsfMaterialEffectsUpdator::updateState(
// Determine the pathCorrection if the material properties exist
pathCorrection =
materialProperties
- ? 1. / fabs(surface->normal().dot(trackParameters->momentum().unit()))
+ ? 1. / std::abs(surface->normal().dot(trackParameters->momentum().unit()))
: 0.;
}
}
@@ -216,7 +216,7 @@ Trk::GsfMaterialEffectsUpdator::preUpdateState(
// Determine the pathCorrection if the material properties exist
pathCorrection =
materialProperties
- ? 1. / fabs(surface->normal().dot(trackParameters->momentum().unit()))
+ ? 1. / std::abs(surface->normal().dot(trackParameters->momentum().unit()))
: 0.;
}
}
@@ -316,7 +316,7 @@ Trk::GsfMaterialEffectsUpdator::postUpdateState(
// Determine the pathCorrection if the material properties exist
pathCorrection =
materialProperties
- ? 1. / fabs(surface->normal().dot(trackParameters->momentum().unit()))
+ ? 1. / std::abs(surface->normal().dot(trackParameters->momentum().unit()))
: 0.;
}
}
@@ -446,7 +446,7 @@ bool
Trk::GsfMaterialEffectsUpdator::updateP(AmgVector(5) & stateVector,
double deltaP) const
{
- double p = 1. / fabs(stateVector[Trk::qOverP]);
+ double p = 1. / std::abs(stateVector[Trk::qOverP]);
p += deltaP;
if (p <= 0.) {
return false;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolution.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolution.cxx
index 438bb317b845b090869ba45e02fbd2a9f18af03d..eabd4ae974bf1d4644bf7deab9d3d72928749d5b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolution.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolution.cxx
@@ -58,6 +58,7 @@ Trk::GsfMaterialMixtureConvolution::finalize()
Trk::MultiComponentState
Trk::GsfMaterialMixtureConvolution::update(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const Trk::MultiComponentState& multiComponentState,
const Trk::Layer& layer,
Trk::PropDirection direction,
@@ -113,6 +114,7 @@ Trk::GsfMaterialMixtureConvolution::update(
QuickCloseComponentsMultiStateMerger::merge(
std::move(cache.multiComponentState), m_maximumNumberOfComponents);
+ ATH_MSG_DEBUG("UPDATE update N in: " << multiComponentState.size() <<" N out: "<< mergedState.size() );
if (mergedState.empty()) {
return {};
}
@@ -128,6 +130,7 @@ Trk::GsfMaterialMixtureConvolution::update(
Trk::MultiComponentState
Trk::GsfMaterialMixtureConvolution::preUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const Trk::MultiComponentState& multiComponentState,
const Trk::Layer& layer,
Trk::PropDirection direction,
@@ -159,7 +162,7 @@ Trk::GsfMaterialMixtureConvolution::preUpdate(
Trk::MultiComponentState updatedState = m_updator->preUpdateState(
*component, layer, direction, particleHypothesis);
-
+ ATH_MSG_DEBUG("PREUPDATE update component result size: " << updatedState.size() );
if (updatedState.empty()) {
continue;
}
@@ -172,11 +175,12 @@ Trk::GsfMaterialMixtureConvolution::preUpdate(
"Component could not be added to the state in the assembler");
}
}
-
+ ATH_MSG_DEBUG("PREUPDATE before merge N: " << cache.multiComponentState.size() );
Trk::MultiComponentState mergedState =
QuickCloseComponentsMultiStateMerger::merge(
std::move(cache.multiComponentState), m_maximumNumberOfComponents);
+ ATH_MSG_DEBUG("PREUPDATE update N in: " << multiComponentState.size() <<" N out: "<< mergedState.size() );
if (mergedState.empty()) {
return {};
}
@@ -192,6 +196,7 @@ Trk::GsfMaterialMixtureConvolution::preUpdate(
Trk::MultiComponentState
Trk::GsfMaterialMixtureConvolution::postUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>&,
const Trk::MultiComponentState& multiComponentState,
const Trk::Layer& layer,
Trk::PropDirection direction,
@@ -243,6 +248,7 @@ Trk::GsfMaterialMixtureConvolution::postUpdate(
QuickCloseComponentsMultiStateMerger::merge(
std::move(cache.multiComponentState), m_maximumNumberOfComponents);
+ ATH_MSG_DEBUG("POSTUPDATE update N in: " << multiComponentState.size() <<" N out: "<< mergedState.size() );
if (mergedState.empty()) {
return {};
}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolutionLM.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolutionLM.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3f8405742c5811a5d415ed0fcc7d8097ba3fd43f
--- /dev/null
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMaterialMixtureConvolutionLM.cxx
@@ -0,0 +1,503 @@
+/*
+ Copyright (C) 2020-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+/*************************************************************************************
+ GsfMaterialMixtureConvolutionLM.cxx - description
+ -------------------------------------------------
+author : amorley
+email : amorley@cern.ch
+decription : Implementation code for GSF material mixture convolution that uses less mem
+************************************************************************************/
+
+#include "TrkGaussianSumFilter/GsfMaterialMixtureConvolutionLM.h"
+#include "TrkGaussianSumFilter/IMultiStateMaterialEffectsUpdator.h"
+#include "TrkGaussianSumFilter/IMultiStateMaterialEffects.h"
+#include "TrkGaussianSumFilter/MultiComponentStateAssembler.h"
+#include "TrkGaussianSumFilter/MultiComponentStateCombiner.h"
+#include "TrkGaussianSumFilter/QuickCloseComponentsMultiStateMerger.h"
+#include "TrkGeometry/Layer.h"
+#include "TrkGeometry/MaterialProperties.h"
+#include "TrkGaussianSumFilter/AlignedDynArray.h"
+#include "TrkGaussianSumFilter/KLGaussianMixtureReduction.h"
+
+#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
+#include "TrkSurfaces/PerigeeSurface.h"
+
+Trk::GsfMaterialMixtureConvolutionLM::GsfMaterialMixtureConvolutionLM(
+ const std::string& type,
+ const std::string& name,
+ const IInterface* parent)
+ : AthAlgTool(type, name, parent)
+{
+ declareInterface<IMaterialMixtureConvolution>(this);
+}
+
+Trk::GsfMaterialMixtureConvolutionLM::~GsfMaterialMixtureConvolutionLM() = default;
+
+StatusCode
+Trk::GsfMaterialMixtureConvolutionLM::initialize()
+{
+
+ ATH_CHECK( m_updator.retrieve() );
+
+ ATH_CHECK( m_materialEffects.retrieve() );
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+Trk::GsfMaterialMixtureConvolutionLM::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
+/* ==========================================
+ Update with full material effects
+ ========================================== */
+
+Trk::MultiComponentState
+Trk::GsfMaterialMixtureConvolutionLM::update(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>& caches,
+ const Trk::MultiComponentState& multiComponentState,
+ const Trk::Layer& layer,
+ Trk::PropDirection direction,
+ Trk::ParticleHypothesis particleHypothesis) const
+{
+
+ const Trk::MaterialProperties* materialProperties =
+ layer.fullUpdateMaterialProperties(*(multiComponentState.begin()->first));
+
+ if (!materialProperties) {
+ ATH_MSG_DEBUG("UPDATE but no material properties!!!");
+ }
+
+ Trk::MultiComponentState updatedMergedState = update( caches,
+ multiComponentState,
+ layer,
+ direction,
+ particleHypothesis,
+ Normal );
+ ATH_MSG_DEBUG("UPDATE update N in: " << multiComponentState.size() <<" N out: "<< updatedMergedState.size() );
+ if (updatedMergedState.empty()) {
+ return {};
+ }
+ // Renormalise state
+ MultiComponentStateHelpers::renormaliseState(updatedMergedState);
+
+ return updatedMergedState;
+}
+
+/* ==========================================
+ Update with pre-update material effects
+========================================== */
+
+Trk::MultiComponentState
+Trk::GsfMaterialMixtureConvolutionLM::preUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>& caches,
+ const Trk::MultiComponentState& multiComponentState,
+ const Trk::Layer& layer,
+ Trk::PropDirection direction,
+ Trk::ParticleHypothesis particleHypothesis) const
+{
+ const Trk::MaterialProperties* materialProperties =
+ layer.fullUpdateMaterialProperties(*(multiComponentState.begin()->first));
+ if (!materialProperties) {
+ ATH_MSG_DEBUG("PREUPDATE but no material properties!!!");
+ }
+ /* -------------------------------------
+ Preliminary checks
+ ------------------------------------- */
+ Trk::MultiComponentState updatedMergedState = update( caches,
+ multiComponentState,
+ layer,
+ direction,
+ particleHypothesis,
+ Preupdate );
+ ATH_MSG_DEBUG("PREUPDATE update N in: " << multiComponentState.size() <<" N out: "<< updatedMergedState.size() );
+ if (updatedMergedState.empty()) {
+ return {};
+ }
+ // Renormalise state
+ MultiComponentStateHelpers::renormaliseState(updatedMergedState);
+
+ return updatedMergedState;
+}
+
+/* ==========================================
+ Update with post-update material effects
+ ========================================== */
+
+Trk::MultiComponentState
+Trk::GsfMaterialMixtureConvolutionLM::postUpdate(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>& caches,
+ const Trk::MultiComponentState& multiComponentState,
+ const Trk::Layer& layer,
+ Trk::PropDirection direction,
+ Trk::ParticleHypothesis particleHypothesis) const
+{
+
+ const Trk::MaterialProperties* materialProperties =
+ layer.fullUpdateMaterialProperties(*(multiComponentState.begin()->first));
+ if (!materialProperties) {
+ ATH_MSG_DEBUG("POSTUPDATE but no material properties!!!");
+ }
+ /* -------------------------------------
+ Preliminary checks
+ ------------------------------------- */
+
+ Trk::MultiComponentState updatedMergedState = update( caches,
+ multiComponentState,
+ layer,
+ direction,
+ particleHypothesis,
+ Postupdate );
+
+
+ ATH_MSG_DEBUG("POSTUPDATE update N in: " << multiComponentState.size() <<" N out: "<< updatedMergedState.size() );
+ if (updatedMergedState.empty()) {
+ return {};
+ }
+ // Renormalise state
+ MultiComponentStateHelpers::renormaliseState(updatedMergedState);
+
+
+ return updatedMergedState;
+}
+
+
+
+Trk::MultiComponentState Trk::GsfMaterialMixtureConvolutionLM::update(
+ std::vector<Trk::IMultiStateMaterialEffects::Cache>& caches,
+ const Trk::MultiComponentState& inputState,
+ const Trk::Layer& layer,
+ Trk::PropDirection direction,
+ Trk::ParticleHypothesis particleHypothesis,
+ MaterialUpdateType updateType) const
+{
+
+
+ // Check the multi-component state is populated
+ if (inputState.empty()) {
+ ATH_MSG_DEBUG("Multi component state passed to update is "
+ "not populated... returning 0");
+ return {};
+ }
+
+
+ double updateFactor(1.);
+ // Full method does this for each component which i don't think this is needed
+ if( updateType == Preupdate ){
+ updateFactor = layer.preUpdateMaterialFactor(*inputState.front().first, direction);
+ ATH_MSG_DEBUG("Material effects update prior to propagation using layer "
+ "information and particle hypothesis: "
+ << particleHypothesis);
+ } else if ( updateType == Postupdate) {
+ updateFactor = layer.postUpdateMaterialFactor( *inputState.front().first, direction);
+ ATH_MSG_DEBUG("Material effects update after propagation using layer "
+ "information and particle hypothesis: "
+ << particleHypothesis);
+ }
+
+ if( updateFactor < 0.01 ){
+ //Bail out as factor is too small to bother about
+ return {};
+ }
+
+
+ caches.resize( inputState.size() );
+
+ // Fill cache and work out how many final components there should be
+ size_t n(0);
+ for( size_t i(0) ; i < inputState.size(); ++i){
+ const AmgSymMatrix(5)* measuredCov = inputState[i].first->covariance();
+ //If the momentum is too dont apply material effects
+ if( inputState[i].first->momentum().mag() <= m_momentumCut ){
+ ATH_MSG_DEBUG("Ignoring material effects... Momentum too low");
+ caches[i].resetAndAddDummyValues();
+ caches[i].deltaParameters[0] = inputState[i].first->parameters();
+ caches[i].weights[0] = inputState[i].second;
+ if (measuredCov) {
+ caches[i].deltaCovariances[0] += *measuredCov;
+ }
+ n += caches[i].weights.size();
+ continue;
+ }
+
+ // Get the material effects and store them in the cache
+ std::pair< const Trk::MaterialProperties*, double > matPropPair =
+ getMaterialProperties( inputState[i].first.get(), layer );
+
+ if( !matPropPair.first ) {
+ ATH_MSG_DEBUG("No material properties .. dont apply material effects");
+ caches[i].resetAndAddDummyValues();
+ caches[i].deltaParameters[0] = inputState[i].first->parameters();
+ caches[i].weights[0] = inputState[i].second;
+ if (measuredCov) {
+ caches[i].deltaCovariances[0] += *measuredCov;
+ }
+ n += caches[i].weights.size();
+ continue;
+ }
+
+ // Apply the update factor
+ matPropPair.second *= updateFactor;
+
+ m_materialEffects->compute(caches[i],
+ inputState[i],
+ *matPropPair.first,
+ matPropPair.second,
+ direction,
+ particleHypothesis);
+
+ // check all vectors have the same size
+ if (caches[i].weights.size() != caches[i].deltaPs.size()) {
+ ATH_MSG_ERROR("Inconsistent number of components in the updator!!! no material effect will be applied");
+ caches[i].resetAndAddDummyValues();
+ }
+
+ for( size_t j(0); j < caches[i].weights.size(); ++j){
+ if (measuredCov) {
+ caches[i].deltaCovariances[j] += *measuredCov;
+ } else {
+ caches[i].deltaCovariances[j].setZero();
+ }
+ caches[i].deltaParameters[j] = inputState[i].first->parameters();
+ // Adjust the momentum of the component's parameters vector here. Check to
+ // make sure update is good.
+ if (!updateP(caches[i].deltaParameters[j][Trk::qOverP], caches[i].deltaPs[j])) {
+ ATH_MSG_ERROR("Cannot update state vector momentum!!! return nullptr");
+ return {};
+ }
+ // Store component weight
+ caches[i].weights[j] *= inputState[i].second;
+ }
+
+ n += caches[i].weights.size();
+ }
+
+
+ // Fill information for to calculate which components to merge
+ // Inaddition scan all components for covariance matrices. If one or more component
+ // is missing an error matrix, component reduction is impossible.
+ bool componentWithoutMeasurement = false;
+
+ GSFUtils::AlignedDynArray<GSFUtils::Component1D, GSFUtils::alignment> components(n);
+ size_t k(0);
+ std::vector< std::pair<size_t,size_t> > indices(n);
+ for( size_t i(0) ; i < inputState.size(); ++i){
+ for( size_t j(0); j < caches[i].weights.size(); ++j ){
+ const AmgSymMatrix(5)* measuredCov = inputState[i].first->covariance();
+ // Fill in infomation
+ const double cov =
+ measuredCov ? caches[i].deltaCovariances[j](Trk::qOverP, Trk::qOverP) : -1.;
+ if( !measuredCov ){
+ componentWithoutMeasurement = true;
+ }
+ components[k].mean = caches[i].deltaParameters[j][Trk::qOverP];
+ components[k].cov = cov;
+ components[k].invCov = cov > 0 ? 1. / cov : 1e10;
+ components[k].weight = caches[i].weights[j] ;
+ indices[k] = {i,j};
+ ++k;
+ }
+ }
+
+
+ if(componentWithoutMeasurement){
+ auto result = std::max_element(components.begin(),
+ components.end(),
+ [](const auto& a, const auto& b){
+ return a.weight < b.weight;
+ });
+ auto index = std::distance(components.begin(), result);
+
+ // Build the first TP
+ size_t stateIndex = indices[index].first;
+ size_t materialIndex = indices[index].second;
+
+ AmgVector(5)& updatedStateVector = caches[stateIndex].deltaParameters[materialIndex];
+ const AmgSymMatrix(5)* measuredCov = inputState[stateIndex].first->covariance();
+ AmgSymMatrix(5)* updatedCovariance = nullptr;
+ if (measuredCov && caches[stateIndex].deltaCovariances.size() > materialIndex) {
+ updatedCovariance = new AmgSymMatrix(5)(
+ caches[stateIndex].deltaCovariances[materialIndex] );
+ }
+ Trk::TrackParameters* updatedTrackParameters =
+ inputState[stateIndex].first->associatedSurface().createTrackParameters(
+ updatedStateVector[Trk::loc1],
+ updatedStateVector[Trk::loc2],
+ updatedStateVector[Trk::phi],
+ updatedStateVector[Trk::theta],
+ updatedStateVector[Trk::qOverP],
+ updatedCovariance);
+
+ Trk::ComponentParameters dummyCompParams( updatedTrackParameters, 1.);
+ Trk::MultiComponentState returnMultiState;
+ returnMultiState.push_back(std::move(dummyCompParams));
+ return returnMultiState;
+ }
+
+
+ // Gather the merges -- order is important -- RHS is smaller than LHS
+ std::vector<std::pair<int32_t, int32_t>> merges;
+ if (n > m_maximumNumberOfComponents)
+ merges = findMerges(components.buffer(), n, m_maximumNumberOfComponents);
+
+ //Merge components
+ MultiComponentStateAssembler::Cache assemblerCache;
+ int nMerges(0);
+ std::vector<bool> isMerged( n, false);
+ for (const auto& mergePair : merges) {
+ const int32_t mini = mergePair.first;
+ const int32_t minj = mergePair.second;
+
+
+
+ // Build the first TP
+ size_t stateIndex = indices[mini].first;
+ size_t materialIndex = indices[mini].second;
+ AmgVector(5)& stateVector = caches[stateIndex].deltaParameters[materialIndex];
+ AmgSymMatrix(5)& measuredCov = caches[stateIndex].deltaCovariances[materialIndex];
+
+ // Build the second TP
+ size_t stateIndex2 = indices[minj].first;
+ size_t materialIndex2 = indices[minj].second;
+
+ ++nMerges;
+
+ const AmgVector(5) firstParameters = stateVector;
+ const double firstWeight = caches[stateIndex].weights[materialIndex];
+
+ Trk::MultiComponentStateCombiner::combineParametersWithWeight( caches[stateIndex].deltaParameters[materialIndex],
+ caches[stateIndex].weights[materialIndex],
+ caches[stateIndex2].deltaParameters[materialIndex2],
+ caches[stateIndex2].weights[materialIndex2] );
+
+ Trk::MultiComponentStateCombiner::combineCovWithWeight( firstParameters,
+ measuredCov,
+ firstWeight,
+ caches[stateIndex2].deltaParameters[materialIndex2],
+ caches[stateIndex2].deltaCovariances[materialIndex2],
+ caches[stateIndex2].weights[materialIndex2] );
+
+ isMerged[minj] = true;
+ caches[stateIndex2].deltaParameters[materialIndex2].setZero();
+ caches[stateIndex2].deltaCovariances[materialIndex2].setZero();
+
+ }
+
+ for( size_t i(0); i < n; ++i ){
+ if( isMerged[i] )
+ continue;
+
+ // Build the TP
+ size_t stateIndex = indices[i].first;
+ size_t materialIndex = indices[i].second;
+ AmgVector(5)& stateVector = caches[stateIndex].deltaParameters[materialIndex];
+ AmgSymMatrix(5)& measuredCov = caches[stateIndex].deltaCovariances[materialIndex];
+
+ Trk::TrackParameters* updatedTrackParameters =
+ inputState[stateIndex].first->associatedSurface().createTrackParameters(
+ stateVector[Trk::loc1],
+ stateVector[Trk::loc2],
+ stateVector[Trk::phi],
+ stateVector[Trk::theta],
+ stateVector[Trk::qOverP],
+ new AmgSymMatrix(5)( measuredCov ) );
+
+ double updatedWeight = caches[stateIndex].weights[materialIndex];
+
+ assemblerCache.multiComponentState.emplace_back( updatedTrackParameters, updatedWeight);
+ assemblerCache.validWeightSum += updatedWeight;
+ }
+
+ if( nMerges + assemblerCache.multiComponentState.size() != n ){
+ ATH_MSG_ERROR("Combining complete but merger size is incompatible: " << n << " " << nMerges << " " << assemblerCache.multiComponentState.size() );
+ }
+
+ // Check all weights
+ Trk::MultiComponentState mergedState =
+ MultiComponentStateAssembler::assembledState(assemblerCache);
+
+ if(mergedState.size() > m_maximumNumberOfComponents)
+ ATH_MSG_ERROR("Merging failed, target size: " << m_maximumNumberOfComponents << " final size: " << mergedState.size());
+ return mergedState;
+}
+
+
+bool Trk::GsfMaterialMixtureConvolutionLM::updateP(double& qOverP, double deltaP) const
+{
+ double p = 1. / std::abs(qOverP);
+ p += deltaP;
+ if (p <= 0.) {
+ return false;
+ }
+ qOverP = qOverP > 0. ? 1. / p : -1. / p;
+ return true;
+}
+
+
+
+std::pair< const Trk::MaterialProperties*, double >
+ Trk::GsfMaterialMixtureConvolutionLM::getMaterialProperties(
+ const Trk::TrackParameters* trackParameters,
+ const Trk::Layer& layer) const
+{
+
+ const Trk::MaterialProperties* materialProperties(nullptr);
+ double pathCorrection(0.);
+
+ // Incorporate the reference material
+
+ if (m_useReferenceMaterial) {
+
+ // Get the surface associated with the parameters
+ const Trk::Surface* surface = &(trackParameters->associatedSurface());
+
+ // Only utilise the reference material if an associated detector element
+ // exists
+ if (surface && surface->associatedDetectorElement()) {
+
+ // Get the layer material properties
+ const Trk::LayerMaterialProperties* layerMaterial =
+ layer.layerMaterialProperties();
+
+ // Assign the material properties
+ materialProperties =
+ layerMaterial ? layerMaterial->fullMaterial(trackParameters->position())
+ : nullptr;
+
+ // Determine the pathCorrection if the material properties exist
+ pathCorrection =
+ materialProperties
+ ? 1. / std::abs(surface->normal().dot(trackParameters->momentum().unit()))
+ : 0.;
+ }
+ }
+
+ // Check that the material properties have been defined - if not define them
+ // from the layer information
+ materialProperties = materialProperties
+ ? materialProperties
+ : layer.fullUpdateMaterialProperties(*trackParameters);
+
+ // Bail out if still no material properties can be found
+ if (!materialProperties) {
+ return {nullptr,0};
+ }
+
+ // Define the path correction
+ pathCorrection =
+ pathCorrection > 0.
+ ? pathCorrection
+ : layer.surfaceRepresentation().pathCorrection(
+ trackParameters->position(), trackParameters->momentum());
+
+
+ // The pathlength ( in mm ) is the path correction * the thickness of the
+ // material
+ double pathLength = pathCorrection * materialProperties->thickness();
+
+ return {materialProperties,pathLength};
+}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMeasurementUpdator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMeasurementUpdator.cxx
index 5300da13fe4504ecf5598e5616d0ecc36bbe5d88..5ffb33af9b5b74fbfde643566eaa7e59c37514df 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMeasurementUpdator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfMeasurementUpdator.cxx
@@ -291,6 +291,8 @@ Trk::GsfMeasurementUpdator::calculateFilterStep(
if (stateBeforeUpdate.empty()) {
ATH_MSG_WARNING("Cannot update multi-state with no components!");
return {};
+ } else {
+ ATH_MSG_DEBUG("calculateFilterStep() starting with : " << stateBeforeUpdate.size() );
}
// Calculate the weight of each component after the measurement
@@ -301,6 +303,8 @@ Trk::GsfMeasurementUpdator::calculateFilterStep(
if (stateWithNewWeights.empty()) {
ATH_MSG_DEBUG("Cacluation of state posterior weights failed... Exiting!");
return {};
+ } else {
+ ATH_MSG_DEBUG("calculateFilterStep() after new weights : " << stateWithNewWeights.size() );
}
// Update each component using the specified updator
@@ -312,9 +316,10 @@ Trk::GsfMeasurementUpdator::calculateFilterStep(
for (; component != stateWithNewWeights.end(); ++component) {
- if (fabs((*component).first->parameters()[Trk::qOverP]) > 0.033333) {
+ if ( stateWithNewWeights.size() > 1 &&
+ std::abs((*component).first->parameters()[Trk::qOverP]) > 0.033333) {
ATH_MSG_DEBUG(
- "About to update component with p<30MeV...skipping component! (2)");
+ "About to update component with p<30 MeV...skipping component! (2)");
continue;
}
Trk::FitQualityOnSurface* componentFitQuality = nullptr;
@@ -381,11 +386,15 @@ Trk::GsfMeasurementUpdator::calculateFilterStep(
}
}
+ ATH_MSG_DEBUG("Assembeler cache size : " << cache.multiComponentState.size() );
+
Trk::MultiComponentState assembledUpdatedState =
MultiComponentStateAssembler::assembledState(cache);
if (assembledUpdatedState.empty()) {
return {};
+ } else {
+ ATH_MSG_DEBUG("Assembeled size : " << assembledUpdatedState.size() );
}
fitQoS = std::make_unique<FitQualityOnSurface>(chiSquared, degreesOfFreedom);
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfSmoother.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfSmoother.cxx
index 5ca92fd3350321ac255d8c5535147013f0dd49c5..4dba395e60173c65ce479704856e22a5ed34b18c 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfSmoother.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfSmoother.cxx
@@ -65,6 +65,7 @@ Trk::GsfSmoother::configureTools(
Trk::SmoothedTrajectory*
Trk::GsfSmoother::fit(const EventContext& ctx,
+ Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
const ForwardTrajectory& forwardTrajectory,
const ParticleHypothesis particleHypothesis,
const Trk::CaloCluster_OnTrack* ccot) const
@@ -249,6 +250,7 @@ Trk::GsfSmoother::fit(const EventContext& ctx,
Trk::MultiComponentState extrapolatedState =
m_extrapolator->extrapolate(ctx,
+ extrapolatorCache,
updatedState,
measurement->associatedSurface(),
Trk::oppositeMomentum,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
index cd299aa50059cbf21d3ca271a9c2a553d4b0c4a9..801690c53a60a8a87205bdceb90076194522e17b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
@@ -14,6 +14,7 @@
#include "TrkGaussianSumFilter/MultiComponentStateModeCalculator.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkSurfaces/Surface.h"
+#include "CxxUtils/phihelper.h"
std::unique_ptr<Trk::TrackParameters>
Trk::MultiComponentStateCombiner::combine(
@@ -41,76 +42,91 @@ Trk::MultiComponentStateCombiner::combineWithWeight(
Trk::ComponentParameters& mergeTo,
const Trk::ComponentParameters& addThis)
{
+ const Trk::TrackParameters* secondParameters = addThis.first.get();
+ combineWithWeight( mergeTo, secondParameters->parameters(), secondParameters->covariance(), addThis.second );
+}
- const Trk::TrackParameters* firstParameters = mergeTo.first.get();
- // Check to see if first track parameters are measured or not
- const AmgSymMatrix(5)* firstMeasuredCov = firstParameters->covariance();
+void
+Trk::MultiComponentStateCombiner::combineWithWeight(
+ Trk::ComponentParameters& mergeTo,
+ const AmgVector(5)& secondParameters, const AmgSymMatrix(5)* secondMeasuredCov, const double secondWeight )
+{
double firstWeight = mergeTo.second;
-
- const Trk::TrackParameters* secondParameters = addThis.first.get();
+ auto trackParameters = mergeTo.first.get();
+ const AmgVector(5)& firstParameters = trackParameters->parameters();
+ AmgVector(5) finalParameters( firstParameters );
// Check to see if first track parameters are measured or not
- const AmgSymMatrix(5)* secondMeasuredCov = secondParameters->covariance();
- double secondWeight = addThis.second;
- double totalWeight = firstWeight + secondWeight;
+ double finalWeight = firstWeight;
+ combineParametersWithWeight( finalParameters, finalWeight,
+ secondParameters, secondWeight );
+
+
+ if( trackParameters->covariance() && secondMeasuredCov )
+ {
+ AmgSymMatrix(5) finalMeasuredCov( *trackParameters->covariance() );
+ combineCovWithWeight( firstParameters, finalMeasuredCov, firstWeight,
+ secondParameters, *secondMeasuredCov, secondWeight );
+
+ mergeTo.first->updateParameters(finalParameters, finalMeasuredCov);
+ mergeTo.second = finalWeight;
+ } else {
+ mergeTo.first->updateParameters(finalParameters, nullptr);
+ mergeTo.second = finalWeight;
+ }
+
+}
+
- AmgVector(5) mean;
- mean.setZero();
- AmgVector(5) parameters = secondParameters->parameters();
+
+void Trk::MultiComponentStateCombiner::combineParametersWithWeight(
+ AmgVector(5)& firstParameters, double& firstWeight,
+ const AmgVector(5)& secondParameters, const double secondWeight )
+{
+
+ double totalWeight = firstWeight + secondWeight;
// Ensure that we don't have any problems with the cyclical nature of phi
// Use first state as reference poin
- double deltaPhi = firstParameters->parameters()[2] - parameters[2];
-
+ double deltaPhi = firstParameters[2] - secondParameters[2];
if (deltaPhi > M_PI) {
- parameters[2] += 2 * M_PI;
+ firstParameters[2] -= 2 * M_PI;
} else if (deltaPhi < -M_PI) {
- parameters[2] -= 2 * M_PI;
+ firstParameters[2] += 2 * M_PI;
}
- mean =
- firstWeight * firstParameters->parameters() + secondWeight * parameters;
- mean /= totalWeight;
+
+ firstParameters =
+ firstWeight * firstParameters + secondWeight * secondParameters;
+ firstParameters /= totalWeight;
// Ensure that phi is between -pi and pi
//
- if (mean[2] > M_PI) {
- mean[2] -= 2 * M_PI;
- } else if (mean[2] < -M_PI) {
- mean[2] += 2 * M_PI;
- }
+ firstParameters[2] = CxxUtils::wrapToPi( firstParameters[2] );
+ firstWeight = totalWeight;
+}
+
+
+
+void Trk::MultiComponentStateCombiner::combineCovWithWeight(
+ const AmgVector(5)& firstParameters, AmgSymMatrix(5)& firstMeasuredCov, const double firstWeight,
+ const AmgVector(5)& secondParameters, const AmgSymMatrix(5)& secondMeasuredCov, const double secondWeight )
+{
+
+ double totalWeight = firstWeight + secondWeight;
// Extract local error matrix: Must make sure track parameters are measured,
// ie have an associated error matrix.
- if (firstMeasuredCov && secondMeasuredCov) {
- AmgSymMatrix(5) covariance;
- AmgSymMatrix(5) covariancePart1;
- covariancePart1.setZero();
- AmgSymMatrix(5) covariancePart2;
- covariancePart2.setZero();
-
- covariancePart1 =
- firstWeight * (*firstMeasuredCov) + secondWeight * (*secondMeasuredCov);
- AmgVector(5) parameterDifference =
- firstParameters->parameters() - parameters;
-
- if (parameterDifference[2] > M_PI) {
- parameterDifference[2] -= 2 * M_PI;
- } else if (parameterDifference[2] < -M_PI) {
- parameterDifference[2] += 2 * M_PI;
- }
- covariancePart2 = firstWeight * secondWeight * parameterDifference *
- parameterDifference.transpose();
- covariance = covariancePart1 / totalWeight +
- covariancePart2 / (totalWeight * totalWeight);
-
- mergeTo.first->updateParameters(mean, covariance);
- mergeTo.second = totalWeight;
- } else {
- mergeTo.first->updateParameters(mean, nullptr);
- mergeTo.second = totalWeight;
- }
+ AmgVector(5) parameterDifference = firstParameters - secondParameters;
+ parameterDifference[2] = CxxUtils::wrapToPi( parameterDifference[2] );
+ parameterDifference /= totalWeight;
+
+ firstMeasuredCov *= firstWeight;
+ firstMeasuredCov += secondWeight * secondMeasuredCov;
+ firstMeasuredCov /= totalWeight;
+ firstMeasuredCov += firstWeight * secondWeight * parameterDifference * parameterDifference.transpose();
+
}
std::unique_ptr<Trk::ComponentParameters>
@@ -227,11 +243,7 @@ Trk::MultiComponentStateCombiner::compute(
// Ensure that phi is between -pi and pi
//
- if (mean[2] > M_PI) {
- mean[2] -= 2 * M_PI;
- } else if (mean[2] < -M_PI) {
- mean[2] += 2 * M_PI;
- }
+ mean[2] = CxxUtils::wrapToPi( mean[2] );
(*covariance) = covariancePart1 / sumW + covariancePart2 / (sumW * sumW);
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
index 6a7df3538a3181e8e1d5a6b71e940b1d7c319e27..ef97513d3e1bd25ee144eb1e35125a9104733437 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
@@ -13,9 +13,10 @@
#include "TrkGaussianSumFilter/MultiComponentStateModeCalculator.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
#include "TrkParameters/TrackParameters.h"
+#include "CxxUtils/phihelper.h"
namespace {
-const double invsqrt2PI = 1. / sqrt(2. * M_PI);
+constexpr double invsqrt2PI = 1. / sqrt(2. * M_PI);
/** bried method to determine the value of the a gaussian distribution at a
* given value */
@@ -165,11 +166,7 @@ Trk::MultiComponentStateModeCalculator::calculateMode(
}
// Ensure that phi is between -pi and pi
if (i == 2) {
- if (modes[i] > M_PI) {
- modes[i] -= 2 * M_PI;
- } else if (modes[i] < -M_PI) {
- modes[i] += 2 * M_PI;
- }
+ modes[i] = CxxUtils::wrapToPi( modes[i] ) ;
}
}
}
@@ -205,9 +202,9 @@ Trk::MultiComponentStateModeCalculator::fillMixture(
// d0=0, z0=1, phi0=2, theta=3, qOverP=4,
double weight = component->second;
double mean = componentParameters->parameters()[parameter[i]];
- // FIXME ATLASRECTS-598 this fabs() should not be necessary... for some
+ // FIXME ATLASRECTS-598 this std::abs() should not be necessary... for some
// reason cov(qOverP,qOverP) can be negative
- double sigma = sqrt(fabs((*measuredCov)(parameter[i], parameter[i])));
+ double sigma = sqrt(std::abs((*measuredCov)(parameter[i], parameter[i])));
// Ensure that we don't have any problems with the cyclical nature of phi
// Use first state as reference point
@@ -254,7 +251,7 @@ Trk::MultiComponentStateModeCalculator::findMode(
double pdfPreviousMode = pdf(previousMode, i, mixture);
if ((pdfMode + pdfPreviousMode) != 0.0) {
- tolerance = fabs(pdfMode - pdfPreviousMode) / (pdfMode + pdfPreviousMode);
+ tolerance = std::abs(pdfMode - pdfPreviousMode) / (pdfMode + pdfPreviousMode);
} else {
return xStart;
}
@@ -288,7 +285,7 @@ Trk::MultiComponentStateModeCalculator::findModeGlobal(
double mode(0);
double maximum(-1);
- double iterate(fabs(mean / 1000));
+ double iterate(std::abs(mean / 1000));
for (double counter(start); counter < end; counter += iterate) {
double value(pdf(counter, i, mixture));
@@ -342,7 +339,7 @@ Trk::MultiComponentStateModeCalculator::findRoot(
e = b - a;
}
- if (fabs(fc) < fabs(fb)) {
+ if (std::abs(fc) < std::abs(fb)) {
ac_equal = true;
a = b;
b = c;
@@ -352,15 +349,15 @@ Trk::MultiComponentStateModeCalculator::findRoot(
fc = fa;
}
- double tol = 0.5 * tolerance * fabs(b);
+ double tol = 0.5 * tolerance * std::abs(b);
double m = 0.5 * (c - b);
- if (fb == 0 || fabs(m) <= tol) {
+ if (fb == 0 || std::abs(m) <= tol) {
result = b;
return true;
}
- if (fabs(e) < tol || fabs(fa) <= fabs(fb)) {
+ if (std::abs(e) < tol || std::abs(fa) <= std::abs(fb)) {
// Bounds decreasing too slowly: use bisection
d = m;
e = m;
@@ -387,8 +384,8 @@ Trk::MultiComponentStateModeCalculator::findRoot(
p = -p;
}
- double min1 = 3 * m * q - fabs(tol * q);
- double min2 = fabs(e * q);
+ double min1 = 3 * m * q - std::abs(tol * q);
+ double min2 = std::abs(e * q);
if (2 * p < (min1 < min2 ? min1 : min2)) {
// Accept the interpolation
e = d;
@@ -403,7 +400,7 @@ Trk::MultiComponentStateModeCalculator::findRoot(
a = b;
fa = fb;
// Evaluate new trial root
- if (fabs(d) > tol) {
+ if (std::abs(d) > tol) {
b += d;
} else {
b += (m > 0 ? +tol : -tol);
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
index 86bd37996bb9a33d5164599b62c9fd34442ab1fe..f53ef90ca49b278f379b7e42a9d7c5a82e9e0e70 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
@@ -1,3 +1,4 @@
+#include "TrkGaussianSumFilter/GsfMaterialMixtureConvolutionLM.h"
#include "TrkGaussianSumFilter/GsfMaterialMixtureConvolution.h"
#include "TrkGaussianSumFilter/GsfCombinedMaterialEffects.h"
#include "TrkGaussianSumFilter/GsfMaterialEffectsUpdator.h"
@@ -8,6 +9,7 @@
#include "TrkGaussianSumFilter/GsfExtrapolator.h"
#include "TrkGaussianSumFilter/GsfSmoother.h"
+DECLARE_COMPONENT( Trk::GsfMaterialMixtureConvolutionLM )
DECLARE_COMPONENT( Trk::GsfMaterialMixtureConvolution )
DECLARE_COMPONENT( Trk::GsfCombinedMaterialEffects )
DECLARE_COMPONENT( Trk::GsfMaterialEffectsUpdator )
diff --git a/Tracking/TrkTruthTracks/TrkTruthTrackTools/src/TruthTrackBuilder.cxx b/Tracking/TrkTruthTracks/TrkTruthTrackTools/src/TruthTrackBuilder.cxx
index d3635339af1c596aaa9e2a97dd37d102f29fee72..b17fa2e045bb4a786da8a7935032b3079cd2d9db 100644
--- a/Tracking/TrkTruthTracks/TrkTruthTrackTools/src/TruthTrackBuilder.cxx
+++ b/Tracking/TrkTruthTracks/TrkTruthTrackTools/src/TruthTrackBuilder.cxx
@@ -237,7 +237,7 @@ Trk::Track* Trk::TruthTrackBuilder::createTrack(const PRD_TruthTrajectory& prdTr
Trk::Track *refittedtrack=m_trackFitter->fit(track,false,materialInteractions);
- //!< @todo : add documentation & find out why we need the fit twice ?
+ //!< Refit a second time to add TRT hits
Trk::Track *refittedtrack2=nullptr;
if (refittedtrack && (int)clusters.size()-i>=9){
Trk::MeasurementSet measset;
@@ -263,7 +263,7 @@ Trk::Track* Trk::TruthTrackBuilder::createTrack(const PRD_TruthTrajectory& prdTr
refittedtrack2=new Trk::Track(refittedtrack->info(),traj2,refittedtrack->fitQuality()->clone());
}
else for (int j=0;j<(int)measset.size();j++) delete measset[j];
- } else {
+ } else if(!refittedtrack){
ATH_MSG_VERBOSE("Track fit of truth trajectory NOT successful, NO track created. ");
return nullptr;
}
@@ -272,7 +272,7 @@ Trk::Track* Trk::TruthTrackBuilder::createTrack(const PRD_TruthTrajectory& prdTr
if (!refittedtrack2 && refittedtrack) refittedtrack2=refittedtrack;
//
- ATH_MSG_WARNING("Track fit of truth trajectory successful, track created. ");
+ ATH_MSG_DEBUG("Track fit of truth trajectory successful, track created. ");
// return what you have
// Before returning, fix the creator
refittedtrack2->info().setPatternRecognitionInfo( Trk::TrackInfo::Pseudotracking);