diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
index 3714b3c87cf03df01420a9f253397973fc210382..fc5c226eb101b73c2dd21c6332ec694aa1aad1c2 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
@@ -12,7 +12,6 @@ decription : Class for fitting according to the Gaussian Sum Filter
formalisation.
********************************************************************************** */
-#include "TrkGaussianSumFilter/IGsfOutlierLogic.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkFitterUtils/FitterTypes.h"
#include "TrkFitterInterfaces/ITrackFitter.h"
@@ -23,9 +22,8 @@ decription : Class for fitting according to the Gaussian Sum Filter
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
-#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/IChronoStatSvc.h"
-#include "xAODEventInfo/EventInfo.h"
#include<atomic>
@@ -36,10 +34,11 @@ class MultiComponentStateOnSurface;
class IMultiComponentStateCombiner;
class IMultiStateExtrapolator;
-class TrackFitInputPreparator;
+class TrackFitInputPreparator;
class IForwardGsfFitter;
class IGsfSmoother;
class Track;
+class FitQuality;
class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
public:
@@ -59,11 +58,11 @@ class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
using ITrackFitter::fit;
/** Refit a track using the Gaussian Sum Filter */
-
+
virtual Track* fit ( const Track&,
const RunOutlierRemoval outlierRemoval = false,
const ParticleHypothesis particleHypothesis = nonInteracting ) const;
-
+
/** Fit a collection of 'PrepRawData' objects using the Gaussian Sum Filter
- This requires that an trackParameters object be supplied also as an initial guess */
virtual Track* fit ( const PrepRawDataSet&,
@@ -84,8 +83,8 @@ class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
const RunOutlierRemoval runOutlier = false,
const ParticleHypothesis matEffects = nonInteracting) const;
- /** Refit a track adding a RIO_OnTrack set
- - This has no form of outlier rejection and will use all hits on orginal track...
+ /** Refit a track adding a RIO_OnTrack set
+ - This has no form of outlier rejection and will use all hits on orginal track...
i.e. very basic impleneation at the moment*/
virtual Track* fit ( const Track&,
const MeasurementSet&,
@@ -106,19 +105,20 @@ class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
/** Produces a perigee from a smoothed trajectory */
const MultiComponentStateOnSurface* makePerigee ( const SmoothedTrajectory*,
const ParticleHypothesis particleHypothesis = nonInteracting ) const;
-
-
+
+ //* Calculate the fit quality */
+ const Trk::FitQuality* buildFitQuality(const Trk::SmoothedTrajectory& ) const;
+
private:
ToolHandle<IMultiStateExtrapolator> m_extrapolator;
-
+
ToolHandle<IMultiStateMeasurementUpdator> m_updator;
ToolHandle<IRIO_OnTrackCreator> m_rioOnTrackCreator;
ToolHandle<IForwardGsfFitter> m_forwardGsfFitter;
ToolHandle<IGsfSmoother> m_gsfSmoother;
- PublicToolHandle<IGsfOutlierLogic> m_outlierLogic
- {this,"GsfOutlierLogic","Trk::GsfOutlierLogic/GsfOutlierLogic",""};
+
bool m_reintegrateOutliers;
@@ -131,12 +131,11 @@ class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
ToolHandle<IMultiComponentStateCombiner> m_stateCombiner;
ServiceHandle<IChronoStatSvc> m_chronoSvc;
-
+
TrackFitInputPreparator* m_inputPreparator;
-
-
- // GSF Fit Statistics
+
+ // GSF Fit Statistics
mutable std::atomic<int> m_FitPRD; // Number of Fit PrepRawData Calls
mutable std::atomic<int> m_FitMeasuremnetBase; // Number of Fit MeasurementBase Calls
mutable std::atomic<int> m_FowardFailure; // Number of Foward Fit Failures
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfOutlierLogic.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfOutlierLogic.h
deleted file mode 100755
index df3bce0820dcbc72516b0ccc9097dda2406966cf..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfOutlierLogic.h
+++ /dev/null
@@ -1,52 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/* *******************************************************************************
- GsfOutlierLogic.h - description
- ---------------------------------
-begin : Tuesday 15th March 2005
-author : atkinson
-email : Tom.Atkinson@cern.ch
-decription : Outlier logic for Gaussian Sum Filter
-********************************************************************************** */
-
-#ifndef TrkGsfOutlierLogic_H
-#define TrkGsfOutlierLogic_H
-
-#include "TrkGaussianSumFilter/IGsfOutlierLogic.h"
-
-#include "TrkFitterUtils/FitterTypes.h"
-
-#include "AthenaBaseComps/AthAlgTool.h"
-
-namespace Trk {
-
-class FitQuality;
-
-class GsfOutlierLogic : public AthAlgTool, virtual public IGsfOutlierLogic {
-
- public:
-
- /** Constructor with AlgTool parameters */
- GsfOutlierLogic( const std::string&, const std::string&, const IInterface* );
-
- /** Virtual destructor */
- virtual ~GsfOutlierLogic() {};
-
- /** AlgTool initialise method */
- StatusCode initialize();
-
- /** AlgTool finalize method */
- StatusCode finalize();
-
- virtual const FitQuality* fitQuality( const SmoothedTrajectory& ) const;
-
- private:
- int m_outputlevel; //!< to cache current output level
-
-};
-
-} // end Trk namespace
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfOutlierLogic.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfOutlierLogic.h
deleted file mode 100755
index fb36bb52e93e9e6ca20009667e00b8b7740ecf75..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfOutlierLogic.h
+++ /dev/null
@@ -1,43 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/* *******************************************************************************
- IGsfOutlierLogic.h - description
- ---------------------------------
-created : Thursday 8th January 2009
-author : amorley
-email : Anthony.Morley@cern.ch
-decription : Abstract interface for outlier logic for Gaussian Sum Filter
-********************************************************************************** */
-
-#ifndef TrkIGsfOutlierLogic_H
-#define TrkIGsfOutlierLogic_H
-
-#include "TrkFitterUtils/FitterTypes.h"
-
-#include "GaudiKernel/IAlgTool.h"
-
-namespace Trk {
-
-class FitQuality;
-
-static const InterfaceID InterfaceID_GsfOutlierLogic ("GsfOutlierLogic", 1, 0);
-
-class IGsfOutlierLogic : virtual public IAlgTool {
-
- public:
-
- /** Virtual destructor */
- virtual ~IGsfOutlierLogic() {};
-
- /** Interface ID method */
- static const InterfaceID& interfaceID() { return InterfaceID_GsfOutlierLogic; };
-
- virtual const FitQuality* fitQuality( const SmoothedTrajectory& ) const = 0;
-
-};
-
-} // end Trk namespace
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateCombiner.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateCombiner.h
index 9bb2af78fee8220d6cc352d604b6365874ce09df..6d401f15d5bd037daa37ca48e1b21650de85b17f 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateCombiner.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateCombiner.h
@@ -8,15 +8,13 @@
created : Thursday 8th January 2009
author : amorley
email : Anthony.Morley@cern.ch
-decription : Abstract interface for the Multi Component State Combiner
+decription : Abstract interface for the Multi Component State Combiner
*******************************************************************************/
#ifndef IMultiComponentStateCombiner_H
#define IMultiComponentStateCombiner_H
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
-#include "TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h"
-
#include "GaudiKernel/IAlgTool.h"
namespace Trk {
@@ -26,21 +24,19 @@ static const InterfaceID IID_MultiComponentStateCombiner("MultiComponentStateCom
class IMultiComponentStateCombiner : virtual public IAlgTool {
public:
-
+
/** Virtual destructor */
virtual ~IMultiComponentStateCombiner () {};
/** AlgTool interface methods */
static const InterfaceID& interfaceID () { return IID_MultiComponentStateCombiner; };
-
+
/** Calculate combined state of many components */
virtual const TrackParameters* combine( const MultiComponentState&, bool useModeTemp = false) const = 0;
/** Calculate combined state and weight of many components */
virtual const ComponentParameters* combineWithWeight( const MultiComponentState&, bool useModeTemp = false ) const = 0;
-
- /** */
- virtual void useMode( bool ) = 0;
+
};
} // end Trk namespace
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h
deleted file mode 100755
index 98a089526724680b753bdda3093f957c016c9ff7..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/***************************************************************************************
- IMultiComponentStateModeCalculator.h - description
- ----------------------------------------------------
-begin : Thursday 6th July 2006
-author : atkinson
-email : Tom.Atkinson@cern.ch
-description : Interface class for AlgTools used to determine the mode (q/p)
- of a gaussian mixture
-***************************************************************************************/
-
-#ifndef Trk_IMultiComponentStateModeCalculator_H
-#define Trk_IMultiComponentStateModeCalculator_H
-
-//#include "GeoPrimitives/GeoPrimitives.h"
-#include "EventPrimitives/EventPrimitives.h"
-#include "GaudiKernel/IAlgTool.h"
-
-namespace Trk{
-
-class MultiComponentState;
-
-//!< Interface ID for IMultiComponentStateModeCalculator
-static const InterfaceID IID_IMultiComponentStateModeCalculator( "IMultiComponentStateModeCalculator", 1, 0 );
-
-class IMultiComponentStateModeCalculator : virtual public IAlgTool{
-
- public:
-
- //!< Virtual destructor
- virtual ~IMultiComponentStateModeCalculator(){};
-
- //!< IAlgTool-AlgTool interface
- static const InterfaceID& interfaceID(){ return IID_IMultiComponentStateModeCalculator; };
-
- //!< Public method to calculate the mode of a multi-component state
- virtual Amg::VectorX calculateMode( const MultiComponentState& ) const = 0;
-
-};
-
-} // end Trk namespace
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
index f0bf8bc07ec99bed54e80c630eebed8de2e78633..9c82c288b73419f4c41c0fa6c39ccb5aa32b567a 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
@@ -19,7 +19,6 @@ description : This class takes a multi-component state and collapses
#include "TrkGaussianSumFilter/IMultiComponentStateCombiner.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
-#include "TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h"
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
@@ -29,7 +28,7 @@ namespace Trk {
class MultiComponentStateCombiner : public AthAlgTool, virtual public IMultiComponentStateCombiner {
public:
-
+
/** Constructor with AlgTool parameters */
MultiComponentStateCombiner (const std::string&, const std::string&, const IInterface*);
@@ -47,26 +46,18 @@ class MultiComponentStateCombiner : public AthAlgTool, virtual public IMultiComp
/** Calculate combined state and weight of many components */
virtual const ComponentParameters* combineWithWeight( const MultiComponentState&, bool useModeTemp = false ) const;
-
- /** */
- virtual void useMode( bool );
-
private:
const ComponentParameters* compute( const MultiComponentState*, bool useModeTemp = false ) const;
-
- ToolHandle<IMultiComponentStateModeCalculator> m_modeCalculator;
-
- mutable bool m_useMode;
- mutable bool m_useModeD0;
- mutable bool m_useModeZ0;
- mutable bool m_useModePhi;
- mutable bool m_useModeTheta;
- mutable bool m_useModeqOverP;
-
- mutable int m_NumberOfCalls;
- mutable float m_fractionPDFused;
+
+ bool m_useMode;
+ bool m_useModeD0;
+ bool m_useModeZ0;
+ bool m_useModePhi;
+ bool m_useModeTheta;
+ bool m_useModeqOverP;
+ float m_fractionPDFused;
};
} // end Trk namespace
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateModeCalculator.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateModeCalculator.h
deleted file mode 100755
index 749ddab0d767c90b749e5e809bdede07664292c4..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateModeCalculator.h
+++ /dev/null
@@ -1,119 +0,0 @@
-/*
- Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
-*/
-
-/***********************************************************************************
- MultiComponentStateModeCalculator.h - description
- ---------------------------------------------------
-begin : Thursday 6th July 2006
-author : atkinson
-email : Tom.Atkinson@cern.ch
-description : Class to calculate the mode (q/p) of a gaussian mixtureArray
-***********************************************************************************/
-
-#ifndef Trk_MultiComponentStateModeCalculator_H
-#define Trk_MultiComponentStateModeCalculator_H
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h"
-#include <atomic>
-#include <array>
-#include <vector>
-
-namespace Trk{
-
-class MultiComponentState;
-
-class MultiComponentStateModeCalculator : public AthAlgTool, virtual public IMultiComponentStateModeCalculator{
-
- public:
-
- //!< Constructor with AlgTool parameters
- MultiComponentStateModeCalculator( const std::string&, const std::string&, const IInterface* );
-
- //!< Virtual destructor
- virtual ~MultiComponentStateModeCalculator();
-
- //!< AlgTool initialise method
- StatusCode initialize();
-
- //!< AlgTool finalise method
- StatusCode finalize();
-
- //!< IMultiComponentStateModeCalculator interface method to calculate mode
- virtual Amg::VectorX calculateMode( const MultiComponentState& ) const override;
-
- private:
-
- //!< Private Mixture structure
- struct Mixture{
-
- // Default constructor
- Mixture():
- weight(0.),
- mean(0.),
- sigma(0.)
- {}
-
- // Constructor with arguments
- Mixture( double aWeight, double aMean, double aSigma ):
- weight( aWeight ),
- mean( aMean ),
- sigma( aSigma )
- {}
- double weight;
- double mean;
- double sigma;
- };
-
-
- //!< Private method to extract the weight, mean and sigma values from the multi-component state
- void fillMixture( std::array<std::vector< Mixture >,5>& mixtureArray,
- const MultiComponentState& ) const;
-
- //!< Private method to find the mode using the Newton-Raphson method based on a starting guess
- double findMode( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double, int ) const;
-
- //!< Private method to determine the pdf of the cashed mixtureArray at a given value
- double pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double, int ) const;
-
- //!< Private method to determine the first order derivative of the pdf at a given value
- double d1pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double, int ) const;
-
- //!< Private method to determine the second order derivative of the pdf at a given value
- double d2pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double, int) const;
-
- //!< Private method to determine the value of the a gaussian distribution at a given value
- double gaus( double x, double mean, double sigma ) const;
-
- double findModeGlobal(const std::array<std::vector< Mixture >,5>& mixtureArray,
- double, int) const;
-
- double width( const std::array<std::vector< Mixture >,5>& mixtureArray,
- int i) const;
-
- double findRoot(const std::array<std::vector< Mixture >,5>& mixtureArray,
- double &result, double xlo, double xhi, double value, double i) const;
-
-
- private:
-
- int m_outputlevel; //!< to cache current output level
-
-
- //ModeCalualtor Stats
- mutable std::atomic<int> m_NumberOfCalls;
- mutable std::atomic<int> m_ConverganceFilures;
- mutable std::atomic<int> m_NoErrorMatrix;
- mutable std::atomic<int> m_MixtureSizeZero;
-
-
-};
-
-} // end Trk namespace
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/CloseComponentsMultiStateMerger.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/CloseComponentsMultiStateMerger.cxx
index f008ec20a67b48c9073dd3976d85f2df096b404c..cdcc903ccf986cd5eab1d0444ccade99d58d1b16 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/CloseComponentsMultiStateMerger.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/CloseComponentsMultiStateMerger.cxx
@@ -27,12 +27,12 @@ Trk::CloseComponentsMultiStateMerger::CloseComponentsMultiStateMerger(
m_stateCombiner("Trk::MultiComponentStateCombiner/CloseComponentsCombiner"),
m_chronoSvc( "ChronoStatSvc", name )
{
-
+
declareInterface<IMultiComponentStateMerger>(this);
declareProperty("MaximumNumberOfComponents", m_maximumNumberOfComponents);
declareProperty("DistanceType", m_distance);
-
+
}
Trk::CloseComponentsMultiStateMerger::~CloseComponentsMultiStateMerger()
@@ -47,7 +47,7 @@ StatusCode Trk::CloseComponentsMultiStateMerger::initialize()
if ( m_chronoSvc.retrieve().isFailure() ) {
msg(MSG::FATAL) << "Failed to retrieve service " << m_chronoSvc << endmsg;
return StatusCode::FAILURE;
- } else
+ } else
msg(MSG::INFO) << "Retrieved service " << m_chronoSvc << endmsg;
// Retrieve the distance tool
@@ -63,8 +63,6 @@ StatusCode Trk::CloseComponentsMultiStateMerger::initialize()
return StatusCode::FAILURE;
}
- m_stateCombiner->useMode(false);
-
// Request an instance of the MultiComponentStateAssembler
if ( m_stateAssembler.retrieve().isFailure() ){
msg(MSG::FATAL) << "Could not retrieve an instance of the mutli-component state assembler... Exiting!" << endmsg;
@@ -92,7 +90,7 @@ const Trk::MultiComponentState* Trk::CloseComponentsMultiStateMerger::merge(cons
// Start the timer
m_chronoSvc->chronoStart("MS::scan");
-
+
if (msgLvl(MSG::VERBOSE)) msg() << "Merging state with " << unmergedState.size() << " components" << endmsg;
bool componentWithoutMeasurement = false;
@@ -166,7 +164,7 @@ const Trk::MultiComponentState* Trk::CloseComponentsMultiStateMerger::merge(cons
while (numberOfComponents > m_maximumNumberOfComponents){
// Reset the merged components multi-map for next iteration
-
+
if ( !mergedComponentsMap.empty() )
mergedComponentsMap.clear();
@@ -186,7 +184,7 @@ const Trk::MultiComponentState* Trk::CloseComponentsMultiStateMerger::merge(cons
// Combine the closest distance components
const Trk::ComponentParameters* combinedComponents = m_stateCombiner->combineWithWeight( *componentsToBeMerged );
-
+
// Erase these components from the unmerged components map. These need to be deleted also as new memory is assigned in the combiner
delete componentsToBeMerged;
@@ -195,7 +193,7 @@ const Trk::MultiComponentState* Trk::CloseComponentsMultiStateMerger::merge(cons
// Insert this merged component into the merged components multi-map
mergedComponentsMap.insert( std::make_pair(combinedComponents->second, *combinedComponents) );
-
+
// Try deleting the minimum distance pair
delete minimumDistancePair;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
index 9dfaa23988fa809c8d57646378187a10f490c0bc..e462befd36e0883b600862d3d65693e1fd57bf92 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
@@ -7,7 +7,7 @@
-------------------------------------
begin : Monday 7th March 2005
author : amorley atkinson
-email : Anthony.Morley@cern.ch Tom.Atkinson@cern.ch
+email : Anthony.Morley@cern.ch Tom.Atkinson@cern.ch
decription : Implementation code for Gaussian Sum Fitter class
********************************************************************************** */
@@ -101,18 +101,15 @@ StatusCode Trk::GaussianSumFitter::initialize()
// Request the GSF smoother - hardwired type and instance name for the GSF
ATH_CHECK( m_gsfSmoother.retrieve() );
- // Request the GSF Outlier m_logic - hardwired type and instance name for the GSF
- ATH_CHECK( m_outlierLogic.retrieve() );
-
// Request the GSF measurement updator - hardwired type and instance name for the GSF
ATH_CHECK ( m_updator.retrieve() );
// Request the GSF extrapolator
ATH_CHECK ( m_extrapolator.retrieve() );
-
+
// Request the state combiner
ATH_CHECK ( m_stateCombiner.retrieve() );
-
+
// Request the RIO_OnTrack creator
// No need to return if RioOnTrack creator tool, only if PrepRawData is used in fit
if( m_rioOnTrackCreator.retrieve().isFailure() ){
@@ -165,7 +162,6 @@ StatusCode Trk::GaussianSumFitter::initialize()
m_inputPreparator = new TrackFitInputPreparator();
-
msg(MSG::INFO) << "Initialisation of " << name() << " was successful" << endmsg;
return StatusCode::SUCCESS;
@@ -183,7 +179,7 @@ StatusCode Trk::GaussianSumFitter::finalize()
msg(MSG::INFO) << "-----------------------------------------------" << endmsg;
msg(MSG::INFO) << " Some Brief GSF Statistics " << endmsg;
msg(MSG::INFO) << "-----------------------------------------------" << endmsg;
-
+
msg(MSG::INFO) << "Number of Fit PrepRawData Calls: "<< m_FitPRD << endmsg;
msg(MSG::INFO) << "Number of Fit MeasurementBase Calls: "<< m_FitMeasuremnetBase << endmsg;
msg(MSG::INFO) << "Number of Forward Fit Failures: "<< m_FowardFailure << endmsg;
@@ -191,7 +187,7 @@ StatusCode Trk::GaussianSumFitter::finalize()
msg(MSG::INFO) << "Number of MakePerigee Failures: "<< m_PerigeeFailure << endmsg;
msg(MSG::INFO) << "Number of Trks that fail fitquality test: "<< m_fitQualityFailure << endmsg;
msg(MSG::INFO) << "-----------------------------------------------" << endmsg;
-
+
msg(MSG::INFO) << "Finalisation of " << name() << " was successful" << endmsg;
return StatusCode::SUCCESS;
@@ -201,7 +197,7 @@ StatusCode Trk::GaussianSumFitter::finalize()
#if 0
StatusCode Trk::GaussianSumFitter::configureTools(const IMultiStateMeasurementUpdator* measurementUpdator, const IRIO_OnTrackCreator* rioOnTrackCreator)
{
-
+
msg(MSG::INFO) << "Configuring the GaussianSumFilter!" << endmsg;
StatusCode sc;
@@ -244,7 +240,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::Track& inputTra
const Trk::RunOutlierRemoval outlierRemoval,
const Trk::ParticleHypothesis particleHypothesis ) const
{
-
+
if (msgLvl(MSG::VERBOSE))
msg() << "Trk::GaussianSumFilter::fit() - Refitting a track" << endmsg;
@@ -278,25 +274,25 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::Track& inputTra
DataVector<const Trk::TrackStateOnSurface>::const_iterator trackStateOnSurface = inputTrack.trackStateOnSurfaces()->begin();
for ( ; trackStateOnSurface != inputTrack.trackStateOnSurfaces()->end(); ++trackStateOnSurface ) {
-
+
if ( !(*trackStateOnSurface) ){
msg(MSG::WARNING) << "This track contains an empty MeasurementBase object that won't be included in the fit" << endmsg;
continue;
}
-
+
if ( (*trackStateOnSurface)->measurementOnTrack() ){
if ( (*trackStateOnSurface)->type(TrackStateOnSurface::Measurement) ){
measurementSet.push_back( (*trackStateOnSurface)->measurementOnTrack() );
}
else if ( m_reintegrateOutliers && (*trackStateOnSurface)->type(TrackStateOnSurface::Outlier) ){
measurementSet.push_back( (*trackStateOnSurface)->measurementOnTrack() );
- }
+ }
}
}
// Apply GSF fit to MeasurementBase objects
return fit( measurementSet, *parametersNearestReference, outlierRemoval, particleHypothesis );
-
+
}
// If refitting of the track is at the PrepRawData level then extract the PrepRawData objects from the input track
@@ -315,19 +311,19 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::Track& inputTra
// Dynamic cast to a RIO_OnTrack object
const Trk::RIO_OnTrack* rioOnTrack = dynamic_cast< const Trk::RIO_OnTrack* >(*measurementOnTrack);
-
+
if ( !rioOnTrack){
msg(MSG::DEBUG) << "Measurement could not be cast as a RIO_OnTrack object... continuing" << endmsg;
continue;
}
-
+
const PrepRawData* prepRawData = rioOnTrack->prepRawData();
-
+
if ( !prepRawData ){
msg(MSG::DEBUG) << "Defined RIO_OnTrack object has no associated PrepRawData object... this object will be ignored in fit" << endmsg;
continue;
}
-
+
prepRawDataSet.push_back( prepRawData );
}
@@ -358,9 +354,8 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
msg() << "Material effects switch: " << particleHypothesis << endmsg;
msg() << "Outlier removal switch: " << outlierRemoval << endmsg;
}
-
++m_FitPRD;
-
+
// Start the timer
Chrono chrono( &(*m_chronoSvc), name() );
@@ -370,7 +365,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
msg(MSG::FATAL) << "PrepRawData set for fit is empty... Exiting!" << endmsg;
return 0;
}
-
+
// A const stl container cannot be sorted. This will re-cast it so that it can.
Trk::PrepRawDataSet sortedPrepRawDataSet = PrepRawDataSet( prepRawDataSet );
@@ -384,7 +379,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
}
- // Perform GSF forwards fit
+ // Perform GSF forwards fit
const ForwardTrajectory* forwardTrajectory = m_forwardGsfFitter->fitPRD( sortedPrepRawDataSet, estimatedParametersNearOrigin, particleHypothesis );
if ( !forwardTrajectory ) {
@@ -392,7 +387,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
++m_FowardFailure;
return 0;
}
-
+
if ( forwardTrajectory->empty() ){
if (msgLvl(MSG::DEBUG)) msg() << "No states in forward trajectory... Exiting!" << endmsg;
++m_FowardFailure;
@@ -402,12 +397,9 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
if (msgLvl(MSG::VERBOSE)) msg() << "*** Forward GSF fit passed! ***" << endmsg;
- // Perform GSF smoother operation
+ // Perform GSF smoother operation
SmoothedTrajectory* smoothedTrajectory = m_gsfSmoother->fit( *forwardTrajectory, particleHypothesis );
-
-
-
// Protect against failed smoother fit
if ( !smoothedTrajectory ) {
if (msgLvl(MSG::DEBUG)) msg() << "Smoother GSF fit failed... Exiting!" << endmsg;
@@ -419,7 +411,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
if (msgLvl(MSG::VERBOSE)) msg() << "*** GSF smoother fit passed! ***" << endmsg;
// Outlier m_logic and track finalisation
- const FitQuality* fitQuality = m_outlierLogic->fitQuality( *smoothedTrajectory );
+ const FitQuality* fitQuality = buildFitQuality( *smoothedTrajectory );
if ( !fitQuality ){
if (msgLvl(MSG::DEBUG)) msg() << "Chi squared could not be calculated... Bailing" << endmsg;
@@ -433,13 +425,13 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
if (outlierRemoval && msgLvl(MSG::DEBUG))
msg(MSG::DEBUG) << "Outlier removal not yet implemented for the Gaussian Sum Filter" << endmsg;
-
-
+
+
if ( m_makePerigee ){
const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface = this->makePerigee( smoothedTrajectory, particleHypothesis );
if (msgLvl(MSG::DEBUG)) msg() << "perigeeMultiStateOnSurface :" << perigeeMultiStateOnSurface << endmsg;
if ( perigeeMultiStateOnSurface ) smoothedTrajectory->push_back( perigeeMultiStateOnSurface );
- else {
+ else {
if (msgLvl(MSG::DEBUG)) msg() << "Perigee asked to be created but failed.....Exiting" << endmsg;
++m_PerigeeFailure;
delete smoothedTrajectory;
@@ -447,24 +439,23 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
return 0;
}
}
-
-
+
+
// Delete forward trajectory. New memory was assigned in ForwardGsfFitter.
delete forwardTrajectory;
-
-
+
//Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
-
-
+
+
// Create new track
- Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
- info.setTrackProperties(TrackInfo::BremFit);
+ Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
+ info.setTrackProperties(TrackInfo::BremFit);
info.setTrackProperties(TrackInfo::BremFitSuccessful);
fittedTrack = new Track(info, smoothedTrajectory, fitQuality );
if ( fittedTrack && msgLvl(MSG::VERBOSE)) {
- msg() << "Fitting of a set of PrepRawData objects is successful" << endmsg;
+ msg() << "Fitting of a set of PrepRawData objects is successful" << endmsg;
msg() << "Track fit chi squared... " << fitQuality->chiSquared() << endmsg;
msg() << "Track fit number of degrees of freedom... " << fitQuality->numberDoF() << endmsg;
}
@@ -497,9 +488,9 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
msg() << "Material effects switch: " << particleHypothesis << endmsg;
msg() << "Outlier removal switch: " << outlierRemoval << endmsg;
}
-
+
++m_FitMeasuremnetBase;
-
+
// Protect against empty PrepRawDataSet object
if ( measurementSet.empty() ) {
msg(MSG::FATAL) << "MeasurementSet for fit is empty... Exiting!" << endmsg;
@@ -509,7 +500,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
// Find the CCOT if it exsists
const Trk::CaloCluster_OnTrack* ccot(0);
Trk::MeasurementSet cleanedMeasurementSet;
-
+
MeasurementSet::const_iterator itSet = measurementSet.begin();
MeasurementSet::const_iterator itSetEnd = measurementSet.end();
for ( ; itSet!=itSetEnd; ++itSet) {
@@ -523,8 +514,8 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
ATH_MSG_DEBUG( "The MeasurementBase object is a Trk::CaloCluster_OnTrack" );
}
}
- }
-
+ }
+
// A const stl container cannot be sorted. This will re-cast it so that it can.
@@ -547,19 +538,19 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
if (msgLvl(MSG::DEBUG)) msg() << "Forward GSF fit failed... Exiting!" << endmsg;
++m_FowardFailure;
return 0;
- }
-
+ }
+
if ( forwardTrajectory->empty() ){
if (msgLvl(MSG::DEBUG)) msg() << "No states in forward trajectory... Exiting!" << endmsg;
delete forwardTrajectory;
++m_FowardFailure;
return 0;
}
-
+
if (msgLvl(MSG::VERBOSE))
msg() << "*** Forward GSF fit passed! ***" << endmsg;
-
+
// Perform GSF smoother operation
SmoothedTrajectory* smoothedTrajectory = m_gsfSmoother->fit( *forwardTrajectory, particleHypothesis, ccot );
@@ -570,17 +561,17 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
delete forwardTrajectory;
return 0;
}
-
+
if (msgLvl(MSG::VERBOSE))
msg() << "*** GSF smoother fit passed! ***" << endmsg;
// Outlier m_logic and track finalisation
- const FitQuality* fitQuality = m_outlierLogic->fitQuality( *smoothedTrajectory );
+ const FitQuality* fitQuality = buildFitQuality( *smoothedTrajectory );
if ( !fitQuality ){
- if (msgLvl(MSG::DEBUG))
+ if (msgLvl(MSG::DEBUG))
msg() << "Chi squared could not be calculated... Bailing" << endmsg;
++m_fitQualityFailure;
delete forwardTrajectory;
@@ -594,9 +585,9 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
if ( m_makePerigee ){
const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface = this->makePerigee( smoothedTrajectory, particleHypothesis );
if (msgLvl(MSG::DEBUG)) msg() << "perigeeMultiStateOnSurface :" << perigeeMultiStateOnSurface << endmsg;
-
+
if ( perigeeMultiStateOnSurface ) smoothedTrajectory->push_back( perigeeMultiStateOnSurface );
- else {
+ else {
if (msgLvl(MSG::DEBUG)) msg() << "Perigee asked to be created but failed.....Exiting" << endmsg;
++m_PerigeeFailure;
delete fitQuality;
@@ -610,16 +601,15 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
delete forwardTrajectory;
-
//Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
-
+
// Create new track
- Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
- info.setTrackProperties(TrackInfo::BremFit);
+ Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
+ info.setTrackProperties(TrackInfo::BremFit);
info.setTrackProperties(TrackInfo::BremFitSuccessful);
- Track* fittedTrack = new Track(info, smoothedTrajectory, fitQuality );
+ Track* fittedTrack = new Track(info, smoothedTrajectory, fitQuality );
if ( fittedTrack ) {
msg(MSG::DEBUG) << "Fitting of a set of MeasurementBase objects is successful" << endmsg;
@@ -628,7 +618,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
} else {
msg(MSG::DEBUG) << "Trk::GaussianSumFilter::fit() failed!" << endmsg;
}
-
+
return fittedTrack;
}
@@ -643,7 +633,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& intrk,
msg() << "--> enter GaussianSumFitter::fit(Track,PrdSet,,)" << endmsg;
msg() << " with Track from author = " << intrk.info().dumpInfo() << endmsg;
}
-
+
// protection, if empty PrepRawDataSet
if (addPrdColl.empty()) {
msg(MSG::WARNING) << "client tries to add an empty PrepRawDataSet to the track fit." << endmsg;
@@ -652,10 +642,10 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& intrk,
/* determine the Track Parameter which is the start of the trajectory,
i.e. closest to the reference point */
- if (msgLvl(MSG::VERBOSE)) msg()<< "get track parameters near origin "
+ if (msgLvl(MSG::VERBOSE)) msg()<< "get track parameters near origin "
<< (m_doHitSorting? "via STL sort" : "from 1st state")
<< endmsg;
-
+
const TrackParameters* estimatedStartParameters = m_doHitSorting
? *(std::min_element(intrk.trackParameters()->begin(),
intrk.trackParameters()->end(),
@@ -664,19 +654,19 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& intrk,
// use external preparator class to prepare PRD set for fitter interface
-
+
Amg::Vector3D referencePosition( m_sortingReferencePoint[0],
m_sortingReferencePoint[1],
m_sortingReferencePoint[2] );
-
+
TrackFitInputPreparator* inputPreparator = new TrackFitInputPreparator(referencePosition);
-
- PrepRawDataSet orderedPRDColl =
+
+ PrepRawDataSet orderedPRDColl =
inputPreparator->stripPrepRawData(intrk,addPrdColl,m_doHitSorting,
true /* do not lose outliers! */);
-
+
delete inputPreparator;
-
+
return fit(orderedPRDColl,*estimatedStartParameters,runOutlier,matEffects);
}
@@ -701,7 +691,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& inputTrack,
msg(MSG::FATAL) << "No estimation of track parameters near origin... Exiting!" << endmsg;
return 0;
}
-
+
// Check that the input track has associated MeasurementBase objects
if ( inputTrack.trackStateOnSurfaces()->empty() ) {
msg(MSG::FATAL) << "Attempting to fit track to empty MeasurementBase collection... Exiting!" << endmsg;
@@ -712,14 +702,14 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& inputTrack,
const Trk::TrackParameters* parametersNearestReference = *( std::min_element( inputTrack.trackParameters()->begin(),
inputTrack.trackParameters()->end(),
*m_trkParametersComparisonFunction ) );
-
-
-
+
+
+
MeasurementSet combinedMS = m_inputPreparator->stripMeasurements (inputTrack, measurementSet,
true, false);
// Apply GSF fit to MeasurementBase objects
- return fit( combinedMS, *parametersNearestReference, runOutlier, matEffects );
+ return fit( combinedMS, *parametersNearestReference, runOutlier, matEffects );
}
Trk::Track* Trk::GaussianSumFitter::fit ( const Track& intrk1,
@@ -728,67 +718,67 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Track& intrk1,
const ParticleHypothesis matEffects) const
{
//Not a great implementation but simple... Just add the hits on track
- // protection against not having measurements on the input tracks
- if (!intrk1.trackStateOnSurfaces() || !intrk2.trackStateOnSurfaces() || intrk1.trackStateOnSurfaces()->size() < 2) {
- msg(MSG::WARNING) << "called to refit empty track or track with too little information, reject fit" << endmsg;
- return 0;
- }
-
- if (!intrk1.trackParameters() || intrk1.trackParameters()->empty()) {
- msg(MSG::WARNING) << "input #1 fails to provide track parameters for seeding the GXF, reject fit" << endmsg;
- return 0;
- }
-
- const TrackParameters* minPar = *intrk1.trackParameters()->begin() ;
+ // protection against not having measurements on the input tracks
+ if (!intrk1.trackStateOnSurfaces() || !intrk2.trackStateOnSurfaces() || intrk1.trackStateOnSurfaces()->size() < 2) {
+ msg(MSG::WARNING) << "called to refit empty track or track with too little information, reject fit" << endmsg;
+ return 0;
+ }
+
+ if (!intrk1.trackParameters() || intrk1.trackParameters()->empty()) {
+ msg(MSG::WARNING) << "input #1 fails to provide track parameters for seeding the GXF, reject fit" << endmsg;
+ return 0;
+ }
+
+ const TrackParameters* minPar = *intrk1.trackParameters()->begin() ;
DataVector<const TrackStateOnSurface>::const_iterator itStates = intrk1.trackStateOnSurfaces()->begin();
DataVector<const TrackStateOnSurface>::const_iterator endState = intrk1.trackStateOnSurfaces()->end();
DataVector<const TrackStateOnSurface>::const_iterator itStates2 = intrk2.trackStateOnSurfaces()->begin();
DataVector<const TrackStateOnSurface>::const_iterator endState2 = intrk2.trackStateOnSurfaces()->end();
Trk::MeasurementSet ms;
-
+
for (;itStates!=endState;++itStates){
if ( !((*itStates)->type(Trk::TrackStateOnSurface::Measurement) || (*itStates)->type(Trk::TrackStateOnSurface::Outlier))) continue;
if (dynamic_cast<const Trk::PseudoMeasurementOnTrack*>((*itStates)->measurementOnTrack())) continue;
ms.push_back((*itStates)->measurementOnTrack());
}
-
+
for (;itStates2!=endState2;++itStates2){
-
+
if ( !((*itStates2)->type(Trk::TrackStateOnSurface::Measurement) || (*itStates2)->type(Trk::TrackStateOnSurface::Outlier))) continue;
-
+
if (dynamic_cast<const Trk::PseudoMeasurementOnTrack*>((*itStates2)->measurementOnTrack())) continue;
-
+
ms.push_back((*itStates2)->measurementOnTrack());
}
-
+
return fit(ms,*minPar,runOutlier,matEffects);
-
+
}
-const Trk::MultiComponentStateOnSurface* Trk::GaussianSumFitter::makePerigee (
+const Trk::MultiComponentStateOnSurface* Trk::GaussianSumFitter::makePerigee (
const Trk::SmoothedTrajectory* smoothedTrajectory,
const Trk::ParticleHypothesis particleHypothesis ) const
{
- if (msgLvl(MSG::VERBOSE))
+ if (msgLvl(MSG::VERBOSE))
msg() << "Trk::GaussianSumFilter::makePerigee... starting" << endmsg;
// Propagate track to perigee
const Trk::PerigeeSurface perigeeSurface;
-
+
const Trk::TrackStateOnSurface* stateOnSurfaceNearestOrigin = smoothedTrajectory->back();
-
+
const Trk::MultiComponentStateOnSurface* multiComponentStateOnSurfaceNearestOrigin = dynamic_cast<const Trk::MultiComponentStateOnSurface*>(stateOnSurfaceNearestOrigin);
-
+
const Trk::MultiComponentState* multiComponentState = 0;
-
+
if ( !multiComponentStateOnSurfaceNearestOrigin ){
- if (msgLvl(MSG::VERBOSE))
+ if (msgLvl(MSG::VERBOSE))
msg() << "State nearest perigee is not a multi-component state... Converting" << endmsg;
-
+
Trk::ComponentParameters componentParameters( stateOnSurfaceNearestOrigin->trackParameters(), 1. );
multiComponentState = new Trk::MultiComponentState( componentParameters );
@@ -796,26 +786,26 @@ const Trk::MultiComponentStateOnSurface* Trk::GaussianSumFitter::makePerigee (
else
multiComponentState = multiComponentStateOnSurfaceNearestOrigin->components();
-
+
// Extrapolate to perigee, taking material effects considerations into account
- const Trk::MultiComponentState* stateExtrapolatedToPerigee = m_extrapolator->extrapolate( *multiComponentState,
- perigeeSurface,
- m_directionToPerigee,
+ const Trk::MultiComponentState* stateExtrapolatedToPerigee = m_extrapolator->extrapolate( *multiComponentState,
+ perigeeSurface,
+ m_directionToPerigee,
false,
particleHypothesis );
-
+
if (!stateExtrapolatedToPerigee){
-
+
if (msgLvl(MSG::DEBUG))
msg() << "Track could not be extrapolated to perigee... returning 0" << endmsg;
-
+
return 0;
}
-
+
// Clean-up & pointer reset
if ( !multiComponentStateOnSurfaceNearestOrigin && stateExtrapolatedToPerigee != multiComponentState )
delete multiComponentState;
-
+
multiComponentState = 0;
// Calculate the mode of the q/p distribution
@@ -832,23 +822,23 @@ const Trk::MultiComponentStateOnSurface* Trk::GaussianSumFitter::makePerigee (
pattern.set(Trk::TrackStateOnSurface::Perigee);
if (fabs(combinedPerigee->parameters()[Trk::qOverP])>1e8) { //GC: protection against 0-momentum track .. this check should NEVER be needed.
- // actual cutoff is 0.01eV track
+ // actual cutoff is 0.01eV track
msg(MSG::ERROR) <<"makePerigee() about to return with 0 momentum!! Returning null instead"<<endmsg;
delete combinedPerigee;
return 0;
}
- const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface
- = new MultiComponentStateOnSurface( 0,
- combinedPerigee,
+ const Trk::MultiComponentStateOnSurface* perigeeMultiStateOnSurface
+ = new MultiComponentStateOnSurface( 0,
+ combinedPerigee,
stateExtrapolatedToPerigee,
- 0,
+ 0,
0,
pattern,
modeQoverP );
msg(MSG::DEBUG) << "makePerigee() returning sucessfully!"<<endmsg;
-
+
return perigeeMultiStateOnSurface;
}
@@ -858,5 +848,31 @@ void Trk::GaussianSumFitter::validationAction() const
}
+const Trk::FitQuality* Trk::GaussianSumFitter::buildFitQuality(const Trk::SmoothedTrajectory& smoothedTrajectory) const
+{
+
+ ATH_MSG_VERBOSE( "Gsf fitQuality" );
+
+ double chiSquared = 0.;
+ int numberDoF = -5;
+ // Loop over all TrackStateOnSurface objects in trajectory
+ SmoothedTrajectory::const_iterator stateOnSurface = smoothedTrajectory.begin();
+ for ( ; stateOnSurface != smoothedTrajectory.end(); ++stateOnSurface ) {
+
+ if ( !(*stateOnSurface)->type(TrackStateOnSurface::Measurement) ) continue;
+ if ( (*stateOnSurface)->fitQualityOnSurface() == 0 ) continue;
+
+ chiSquared += (*stateOnSurface)->fitQualityOnSurface()->chiSquared();
+ numberDoF += (*stateOnSurface)->fitQualityOnSurface()->numberDoF();
+
+ }
+
+ if ( std::isnan( chiSquared ) || chiSquared <= 0. ) return 0;
+
+ const FitQuality* fitQuality = new FitQuality( chiSquared, numberDoF );
+
+ return fitQuality;
+
+}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfOutlierLogic.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfOutlierLogic.cxx
deleted file mode 100755
index 737d84072e4916404870b3417d206d1b1530be7e..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfOutlierLogic.cxx
+++ /dev/null
@@ -1,76 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/* *******************************************************************************
- GsfOutlierLogic.cxx - description
- -----------------------------------
-begin : Tuesday 15th March 2005
-author : amorley, atkinson
-email : Anthony.Morley@cern.ch, Tom.Atkinson@cern.ch
-decription : Implementation code for GsfOutlierLogic class
-********************************************************************************** */
-
-#include "TrkGaussianSumFilter/GsfOutlierLogic.h"
-
-#include "TrkMultiComponentStateOnSurface/MultiComponentStateOnSurface.h"
-#include "TrkEventPrimitives/FitQuality.h"
-
-Trk::GsfOutlierLogic::GsfOutlierLogic(const std::string& type, const std::string& name, const IInterface* parent)
- :
- AthAlgTool(type, name, parent),
- m_outputlevel(0)
-{
-
- declareInterface<IGsfOutlierLogic>(this);
-
-}
-
-StatusCode Trk::GsfOutlierLogic::initialize()
-{
-
- m_outputlevel = msg().level()-MSG::DEBUG; // save the threshold for debug printout in private member
-
- msg(MSG::INFO) << "Initialisation of " << name() << " was successful" << endmsg;
-
- return StatusCode::SUCCESS;
-
-}
-
-StatusCode Trk::GsfOutlierLogic::finalize()
-{
-
- msg(MSG::INFO) << "Finalisation of " << name() << " was successful" << endmsg;
-
- return StatusCode::SUCCESS;
-
-}
-
-const Trk::FitQuality* Trk::GsfOutlierLogic::fitQuality(const Trk::SmoothedTrajectory& smoothedTrajectory) const
-{
-
- msg(MSG::VERBOSE) << "GsfOutlierLogic fitQuality" << endmsg;
-
- double chiSquared = 0.;
- int numberDoF = -5;
-
- // Loop over all TrackStateOnSurface objects in trajectory
- SmoothedTrajectory::const_iterator stateOnSurface = smoothedTrajectory.begin();
-
- for ( ; stateOnSurface != smoothedTrajectory.end(); ++stateOnSurface ) {
-
- if ( !(*stateOnSurface)->type(TrackStateOnSurface::Measurement) ) continue;
- if ( (*stateOnSurface)->fitQualityOnSurface() == 0 ) continue;
-
- chiSquared += (*stateOnSurface)->fitQualityOnSurface()->chiSquared();
- numberDoF += (*stateOnSurface)->fitQualityOnSurface()->numberDoF();
-
- }
-
- if ( std::isnan( chiSquared ) || chiSquared <= 0. ) return 0;
-
- const FitQuality* fitQuality = new FitQuality( chiSquared, numberDoF );
-
- return fitQuality;
-
-}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
index 192f128c3889711e2764de012eaf4995d8aff89a..fbcd111b8e0ee298c10963a0cbc31f2b8c058291 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
@@ -8,62 +8,53 @@
begin : Monday 20th December 2004
author : atkinson
email : Tom.Atkinson@cern.ch
-description : Implementation code for MultiComponentStateCombiner class
+description : Implementation code for MultiComponentStateCombiner class
*********************************************************************************/
#include "TrkGaussianSumFilter/MultiComponentStateCombiner.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkSurfaces/Surface.h"
+#include "MultiComponentStateModeCalculator.h"
+
Trk::MultiComponentStateCombiner::MultiComponentStateCombiner (const std::string& type, const std::string& name, const IInterface* parent)
:
AthAlgTool(type, name, parent),
- m_modeCalculator("Trk::MultiComponentStateModeCalculator/MultiComponentStateModeCalculator"),
m_useMode(false),
m_useModeD0(true),
m_useModeZ0(true),
m_useModePhi(true),
m_useModeTheta(true),
m_useModeqOverP(true),
- m_NumberOfCalls(0),
m_fractionPDFused(1.0)
{
declareInterface<IMultiComponentStateCombiner>(this);
- declareProperty("UseMode", m_useMode );
- declareProperty("ModeCalculator", m_modeCalculator);
+ declareProperty("UseMode", m_useMode, "Calculate mode for all mergers (not recommended)");
declareProperty("UseModeqOverP", m_useModeqOverP );
declareProperty("UseModeD0", m_useModeD0 );
declareProperty("UseModeZ0", m_useModeZ0 );
declareProperty("UseModePhi", m_useModePhi );
declareProperty("UseModeTheta", m_useModeTheta );
declareProperty("FractionPDFused" , m_fractionPDFused);
-
+
}
StatusCode Trk::MultiComponentStateCombiner::initialize()
{
- // Request the mode calculator
- if ( m_modeCalculator.retrieve().isFailure() ){
- ATH_MSG_FATAL( "Unable to retrieve the mode calculator... Exiting!" );
- return StatusCode::FAILURE;
- }
-
- m_NumberOfCalls = 0;
-
if (m_fractionPDFused < 0.1 ){
- ATH_MSG_INFO("Fraction of PDF is set too low begin reset to 1");
- m_fractionPDFused = 1;
+ ATH_MSG_INFO("Fraction of PDF is set too low begin reset to 1");
+ m_fractionPDFused = 1;
}
if (m_fractionPDFused > 1 ){
- ATH_MSG_INFO("Fraction of PDF is set high low begin reset to 1");
- m_fractionPDFused = 1;
+ ATH_MSG_INFO("Fraction of PDF is set high low begin reset to 1");
+ m_fractionPDFused = 1;
}
-
- if (msgLvl(MSG::VERBOSE)) ATH_MSG_VERBOSE( "Initialisation of " << name() << " was successful" );
+
+ ATH_MSG_VERBOSE( "Initialisation of " << name() << " was successful" );
return StatusCode::SUCCESS;
@@ -75,15 +66,14 @@ StatusCode Trk::MultiComponentStateCombiner::finalize()
ATH_MSG_INFO("-----------------------------------------------");
ATH_MSG_INFO(" GSF MCS Combiner Statistics ");
ATH_MSG_INFO("-----------------------------------------------");
- ATH_MSG_INFO("Number of Calls " << m_NumberOfCalls );
ATH_MSG_INFO("Finalisation of " << name() << " was successful" );
-
+
return StatusCode::SUCCESS;
}
const Trk::TrackParameters* Trk::MultiComponentStateCombiner::combine( const Trk::MultiComponentState& uncombinedState, bool useModeTemp ) const
-{
+{
const Trk::ComponentParameters* combinedComponent = compute(&uncombinedState, useModeTemp);
const Trk::TrackParameters* trackParameters = combinedComponent->first;
delete combinedComponent;
@@ -101,7 +91,6 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::combineWithWei
const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const Trk::MultiComponentState* uncombinedState, bool useModeTemp ) const
{
- ++m_NumberOfCalls;
if ( uncombinedState->empty() ){
ATH_MSG_WARNING( "Trying to collapse state with zero components" );
return 0;
@@ -111,10 +100,10 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
// Check to see if first track parameters are measured or not
const AmgSymMatrix(5)* firstMeasuredCov = firstParameters->covariance();
-
+
if ( uncombinedState->size() == 1 )
return new Trk::ComponentParameters(uncombinedState->front().first->clone(), uncombinedState->front().second);
-
+
double sumW(0.);
const int dimension = (uncombinedState->front()).first->parameters().rows();
if (dimension!=5){
@@ -138,7 +127,7 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
Trk::MultiComponentState::const_iterator component = uncombinedState->begin();
double totalWeight(0.);
- for ( ; component != uncombinedState->end() ; ++component){
+ for ( ; component != uncombinedState->end() ; ++component){
double weight = (*component).second;
totalWeight += weight;
}
@@ -146,19 +135,19 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
component = uncombinedState->begin();
- for ( ; component != uncombinedState->end() ; ++component){
+ for ( ; component != uncombinedState->end() ; ++component){
const TrackParameters* trackParameters = (*component).first;
double weight = (*component).second;
-
+
AmgVector(5) parameters = trackParameters->parameters();
-
-
+
+
//Ensure that we don't have any problems with the cyclical nature of phi
//Use first state as reference poin
double deltaPhi = (*uncombinedState->begin()).first->parameters()[2] - parameters[2];
-
+
if( deltaPhi > M_PI ){
parameters[2] += 2 * M_PI;
} else if ( deltaPhi < -M_PI ){
@@ -168,15 +157,15 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
sumW += weight;
mean += weight * parameters;
-
+
// Extract local error matrix: Must make sure track parameters are measured, ie have an associated error matrix.
const AmgSymMatrix(5)* measuredCov = trackParameters->covariance();
-
+
if (measuredCov){
// Changed from errorMatrixInMeasurementFrame
-
+
covariancePart1 += weight * (*measuredCov);
-
+
/* ============================================================================
Loop over all remaining components to find the second part of the covariance
============================================================================ */
@@ -184,17 +173,17 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
Trk::MultiComponentState::const_iterator remainingComponentIterator = component;
for ( ; remainingComponentIterator != uncombinedState->end(); ++remainingComponentIterator) {
-
+
if ( remainingComponentIterator == component ) continue;
-
+
AmgVector(5) parameterDifference = parameters - ((*remainingComponentIterator).first)->parameters();
-
+
double remainingComponentIteratorWeight = (*remainingComponentIterator).second;
-
+
AmgSymMatrix(5) unity;
for(int i(0); i<5;++i){
for(int j(0); j<5;++j){
- unity(i,j) = parameterDifference(i) * parameterDifference(j) ;
+ unity(i,j) = parameterDifference(i) * parameterDifference(j) ;
}
}
covariancePart2 += weight * remainingComponentIteratorWeight * unity ;
@@ -202,33 +191,33 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
} // end loop over remaining components
} // end clause if errors are involved
- if( weight/totalWeight > m_fractionPDFused) break;
+ if( weight/totalWeight > m_fractionPDFused) break;
} // end loop over all components
-
+
mean /= sumW;
-
- //Ensure that phi is between -pi and pi
+
+ //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;
- }
+ }
(*covariance) = covariancePart1 / sumW + covariancePart2 / (sumW * sumW);
-
+
if ( m_useMode || useModeTemp ){
if ( dimension == 5 ){
-
+
// Calculate the mode of the q/p distribution
Amg::VectorX modes(10); modes.setZero();
- modes = m_modeCalculator->calculateMode( *uncombinedState );
+ modes = Trk::MultiComponentStateModeCalculator::calculateMode( *uncombinedState, msgStream() );
if ( msgLvl(MSG::VERBOSE) && modes[4] )
ATH_MSG_VERBOSE( "Calculated mode q/p is: " << modes[4] );
-
+
// Replace mean with mode if qOverP mode is not 0
if (modes[4] != 0){
if (m_useModeD0){
@@ -294,7 +283,7 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
}
} else {
- if (msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG( " Dimension != 5 not updating q/p to mode q/p");
+ ATH_MSG_DEBUG( " Dimension != 5 not updating q/p to mode q/p");
}
@@ -322,9 +311,3 @@ const Trk::ComponentParameters* Trk::MultiComponentStateCombiner::compute( const
return combinedComponentParameters;
}
-
-
-void Trk::MultiComponentStateCombiner::useMode( bool useModeTemp )
-{
- m_useMode = useModeTemp;
-}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
index 6feba1f8876e81113bb7223af3d61d6d3f70ff51..c51529a7eaaf7b5d89cbc71ed3944a82be0d3455 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx
@@ -11,146 +11,100 @@ email : amorley@cern.ch
description : Implementation code for MultiComponentStateModeCalculator class
****************************************************************************************/
-#include "TrkGaussianSumFilter/MultiComponentStateModeCalculator.h"
+#include "MultiComponentStateModeCalculator.h"
#include "TrkMultiComponentStateOnSurface/MultiComponentState.h"
-//#include "EventPrimitives/AmgMatrixPlugin.h"
-
#include "TrkParameters/TrackParameters.h"
#include <map>
-Trk::MultiComponentStateModeCalculator::MultiComponentStateModeCalculator( const std::string& type, const std::string& name, const IInterface* parent )
- :
- AthAlgTool( type, name, parent ),
- m_outputlevel(1),
- m_NumberOfCalls(0),
- m_ConverganceFilures(0),
- m_NoErrorMatrix(0),
- m_MixtureSizeZero(0)
-{
+#include "CxxUtils/AthUnlikelyMacros.h"
- declareInterface<IMultiComponentStateModeCalculator>(this);
-
-}
-Trk::MultiComponentStateModeCalculator::~MultiComponentStateModeCalculator()
-{}
-StatusCode Trk::MultiComponentStateModeCalculator::initialize()
+Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode( const Trk::MultiComponentState& multiComponentState, MsgStream &log )
{
- m_outputlevel = msg().level()-MSG::DEBUG; // save the threshold for debug printout in private member
-
- m_NumberOfCalls = 0;
- m_ConverganceFilures = 0;
- m_NoErrorMatrix = 0;
- m_MixtureSizeZero = 0;
-
- msg(MSG::INFO) << "Initialisation of " << name() << " was successful" << endmsg;
-
- return StatusCode::SUCCESS;
-
-}
-
-StatusCode Trk::MultiComponentStateModeCalculator::finalize()
-{
-
- msg(MSG::INFO) << "-----------------------------------------------"<< endmsg;
- msg(MSG::INFO) << " GSF ModeCalulator Statistics "<< endmsg;
- msg(MSG::INFO) << "-----------------------------------------------"<< endmsg;
- msg(MSG::INFO) << "Number of Calls " << m_NumberOfCalls << endmsg;
- msg(MSG::INFO) << "No Error Matrix " << m_NoErrorMatrix << endmsg;
- msg(MSG::INFO) << "MixtureSize = Zero " << m_MixtureSizeZero << endmsg;
- msg(MSG::INFO) << "Failed to converge " << m_ConverganceFilures << endmsg;
- msg(MSG::INFO) << "Finalisation of " << name() << " was successful" << endmsg;
-
- return StatusCode::SUCCESS;
-
-}
-
-Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode(const Trk::MultiComponentState& multiComponentState ) const
-{
- ++m_NumberOfCalls;
-
Amg::VectorX modes(10) ;
modes.setZero();
-
+
// Check to see if the multi-component state is measured
if ( !multiComponentState.isMeasured() ){
- ATH_MSG_DEBUG( "Mixture has no error matricies... Exiting ");
- ++m_NoErrorMatrix;
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) ) log << "Mixture has no error matricies... Exiting " << endmsg;
return modes;
}
- std::array<std::vector< Mixture >,5> mixtureArray{};
- fillMixture( mixtureArray,multiComponentState );
-
- int tmp_MixtureSizeZero{0};
+
+ std::array<std::vector< Mixture >,5> mixture;
+
+ fillMixture( multiComponentState, mixture );
+
for (int i=0 ; i<5 ; i++){
std::multimap< double, double, std::greater<double> > para_startMap;
- // Loop over the mixtureArray and find the starting point
- std::vector<Mixture>::const_iterator component = mixtureArray[i].begin();
-
- if (mixtureArray[i].size() == 0) {
- ++tmp_MixtureSizeZero;
- }
-
- for( ; component != mixtureArray[i].end() ; ++component ){
- para_startMap.insert( std::pair<double, double>( pdf( mixtureArray,component->mean, i), component->mean ) );
- }
-
+
+ // Loop over the mixture and find the starting point
+ auto component = mixture[i].begin();
+
+ for( ; component != mixture[i].end() ; ++component )
+ para_startMap.insert( std::pair<double, double>( pdf( component->mean, i, mixture), component->mean ) );
+
double para_start = para_startMap.begin()->second;
- modes[i] = findMode( mixtureArray,para_start, i ) ;
- // Calculate the FWHM and return this back so that it can be used to correct the covariance matrix
+ modes[i] = findMode( para_start, i, mixture, log ) ;
+ // Calculate the FWHM and return this back so that it can be used to correct the covariance matrix
if( para_start != modes[i] ){
// mode calculation was successful now calulate FWHM
- double currentWidth = width(mixtureArray,i);
+ double currentWidth = width(i, mixture);
modes[i+5] = -1; // Failure is flagged with a value less than 0;
-
- double pdfVal = pdf( mixtureArray,modes[i], i );
- double highX(0);
- double lowX(0);
+
+ double pdfVal = pdf( modes[i], i, mixture );
+ double highX(0);
+ double lowX(0);
double upperbound =modes[i] + 1.5 * currentWidth;
while(true){
- if( pdf( mixtureArray,upperbound, i ) > pdfVal*0.5 ){
- upperbound += currentWidth;
+ if( pdf( upperbound, i, mixture ) > pdfVal*0.5 ){
+ upperbound += currentWidth;
} else {
break;
}
}
- ATH_MSG_VERBOSE ( "HighX PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,modes[i], i ) << ", " << pdf( mixtureArray,upperbound, i ) << "]") ;
+ if( ATH_UNLIKELY( log.level() <= MSG::VERBOSE ) )
+ log << "HighX PDFval, high val, low val [ " << pdfVal << ", " << pdf( modes[i], i, mixture ) << ", " << pdf( upperbound, i, mixture ) << "]" << endmsg;
+
+ bool highXFound = findRoot( highX, modes[i], upperbound, pdfVal*0.5, i, mixture, log );
- bool highXFound = findRoot( mixtureArray,highX, modes[i], upperbound, pdfVal*0.5, i);
-
double lowerbound =modes[i] - 1.5 * currentWidth;
while(true){
- if( pdf( mixtureArray,lowerbound, i ) > pdfVal*0.5 ){
- lowerbound -= currentWidth;
+ if( pdf( lowerbound, i, mixture ) > pdfVal*0.5 ){
+ lowerbound -= currentWidth;
} else {
break;
}
}
- ATH_MSG_VERBOSE ( "LowX PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,lowerbound, i ) << ", " << pdf( mixtureArray,modes[i], i ) << "]") ;
+ if( ATH_UNLIKELY( log.level() <= MSG::VERBOSE ) )
+ log << "LowX PDFval, high val, low val [ " << pdfVal << ", " << pdf( lowerbound, i, mixture ) << ", " << pdf( modes[i], i, mixture ) << "]" << endmsg;
+
+ bool lowXFound = findRoot( lowX , lowerbound, modes[i], pdfVal*0.5, i, mixture, log);
- bool lowXFound = findRoot( mixtureArray,lowX , lowerbound, modes[i], pdfVal*0.5, i);
-
if (highXFound && lowXFound ){
double FWHM = highX - lowX;
- ATH_MSG_DEBUG ( "PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,highX, i ) << ", " << pdf( mixtureArray,lowX, i ) << "]") ;
-
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "PDFval, high val, low val [ " << pdfVal << ", " << pdf( highX, i, mixture ) << ", " << pdf( lowX, i, mixture) << "]" << endmsg;
if( FWHM <= 0 ) {
- ATH_MSG_DEBUG(i << " Width is neagtive? " << highX << " " << lowX << " " << modes[i] );
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << i << " Width is neagtive? " << highX << " " << lowX << " " << modes[i] << endmsg;
} else {
- ATH_MSG_DEBUG(i << " Width is positive? " << highX << " " << lowX << " " << modes[i] );
- ATH_MSG_DEBUG("Old & New width " << currentWidth << " " << FWHM/2.35 << " High side only " << 2*( highX - modes[i])/2.355 );
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) ){
+ log << i << " Width is positive? " << highX << " " << lowX << " " << modes[i] << endmsg;
+ log << "Old & New width " << currentWidth << " " << FWHM/2.35 << " High side only " << 2*( highX - modes[i])/2.355 << endmsg;
+ }
modes[i+5] = FWHM/2.35482; //2 * sqrt( 2* log(2))
}
} else {
- ATH_MSG_DEBUG( i << " Failed to find 1/2 width " );
-
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) ){
+ log << i << " Failed to find 1/2 width " << endmsg;
+ }
}
//Ensure that phi is between -pi and pi
if(i==2){
@@ -163,59 +117,53 @@ Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode(const Trk::Mu
}
para_startMap.clear();
}
- m_MixtureSizeZero+=tmp_MixtureSizeZero;
return modes;
}
-void Trk::MultiComponentStateModeCalculator::fillMixture( std::array<std::vector< Mixture >,5>& mixtureArray,
- const Trk::MultiComponentState& multiComponentState ) const
+void Trk::MultiComponentStateModeCalculator::fillMixture( const Trk::MultiComponentState& multiComponentState, std::array<std::vector< Mixture >,5>& mixture )
{
-
+
for (int i=0; i<5; i++){
- mixtureArray[i].clear();
+ mixture[i].clear();
}
-
+
// Loop over all the components in the multi-component state
Trk::MultiComponentState::const_iterator component = multiComponentState.begin();
Trk::ParamDefs parameter[5]={Trk::d0,Trk::z0,Trk::phi,Trk::theta,Trk::qOverP};
for( ; component != multiComponentState.end() ; ++component ){
- for (int i=0; i<5; ++i){
- const Trk::TrackParameters* componentParameters = component->first;
-
- const AmgSymMatrix(5)* measuredCov = componentParameters->covariance();
-
- if ( !measuredCov ) return;
-//Enums for Perigee //
-// 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 reason cov(qOverP,qOverP) can be negative
- double sigma = sqrt(fabs((*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
- if(i==2){ //phi
- double deltaPhi = multiComponentState.begin()->first->parameters()[2] -mean;
- if( deltaPhi > M_PI ){
- mean += 2 * M_PI;
- } else if ( deltaPhi < -M_PI ){
- mean -= 2 * M_PI;
+ for (int i=0; i<5; ++i){
+ const Trk::TrackParameters* componentParameters = component->first;
+
+ const AmgSymMatrix(5)* measuredCov = componentParameters->covariance();
+
+ if ( !measuredCov ) return;
+ //Enums for Perigee //
+ // 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 reason cov(qOverP,qOverP) can be negative
+ double sigma = sqrt(fabs((*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
+ if(i==2){ //phi
+ double deltaPhi = multiComponentState.begin()->first->parameters()[2] -mean;
+ if( deltaPhi > M_PI ){
+ mean += 2 * M_PI;
+ } else if ( deltaPhi < -M_PI ){
+ mean -= 2 * M_PI;
+ }
}
+ Mixture mix(weight, mean, sigma );
+ mixture[i].push_back( mix );
}
-
- Mixture newmixture(weight, mean, sigma );
- mixtureArray[i].push_back( newmixture);
- }
}
-
return;
-
}
-double Trk::MultiComponentStateModeCalculator::findMode( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double xStart, int i ) const
+double Trk::MultiComponentStateModeCalculator::findMode( double xStart, int i, std::array<std::vector< Mixture >,5>& mixture, MsgStream &log )
{
int iteration(0);
@@ -224,45 +172,48 @@ double Trk::MultiComponentStateModeCalculator::findMode( const std::array<std::v
double mode( xStart );
- while( iteration < 20 && tolerance > 1.e-8 ){
+ while( iteration < 20 && tolerance > 1.e-8 ){
- double previousMode(mode);
-
- if (d2pdf( mixtureArray,previousMode, i) != 0.0) {
- mode = previousMode - d1pdf(mixtureArray,previousMode, i ) / d2pdf( mixtureArray,previousMode, i );
- }
- else {
- ATH_MSG_DEBUG( "Second derivative is zero ... Returning the original value" );
- return xStart;
- }
+ double previousMode(mode);
+ double d2pdfVal = d2pdf( previousMode, i, mixture);
- if ( ( pdf(mixtureArray,mode, i) + pdf(mixtureArray,previousMode, i)) != 0.0 ) {
- tolerance = fabs( pdf(mixtureArray,mode, i) - pdf(mixtureArray,previousMode, i) )
- / ( pdf(mixtureArray,mode, i) + pdf(mixtureArray,previousMode, i) );
- }
- else {
- ATH_MSG_DEBUG( "Dividing by zero ... Returning the original value" );
- return xStart;
- }
+ if (d2pdfVal != 0.0) {
+ mode = previousMode - d1pdf( previousMode, i, mixture ) / d2pdfVal;
+ } else {
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "Second derivative is zero ... Returning the original value" << endmsg;
+ return xStart;
+ }
- ++iteration;
- }
+ double pdfMode = pdf(mode, i, mixture);
+ double pdfPreviousMode = pdf(previousMode, i, mixture);
- if ( iteration >= 20 ){
- ATH_MSG_DEBUG( "Could not converge to the mode within allowed iterations... Returning the original value" );
- ++m_ConverganceFilures;
+ if ( (pdfMode + pdfPreviousMode) != 0.0 ) {
+ tolerance = fabs( pdfMode - pdfPreviousMode ) / ( pdfMode + pdfPreviousMode );
+ } else {
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "Dividing by zero ... Returning the original value" << endmsg;
return xStart;
}
- return mode;
+ ++iteration;
+
+ }
+
+ if ( iteration >= 20 ){
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "Could not converge to the mode within allowed iterations... Returning the original value" << endmsg;
+ return xStart;
+ }
+
+ return mode;
}
-double Trk::MultiComponentStateModeCalculator::findModeGlobal( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double mean,int i) const
+double Trk::MultiComponentStateModeCalculator::findModeGlobal(double mean,int i, std::array<std::vector< Mixture >,5>& mixture)
{
-
+
double start(-1);
double end(1);
if (mean > 0.0) {
@@ -277,9 +228,9 @@ double Trk::MultiComponentStateModeCalculator::findModeGlobal( const std::array<
double mode(0);
double maximum(-1);
double iterate(fabs(mean/1000));
-
+
for (double counter(start); counter < end; counter+=iterate) {
- double value(pdf(mixtureArray,counter,i));
+ double value(pdf(counter,i, mixture));
if (value > maximum) {
maximum = value;
mode = counter;
@@ -288,30 +239,28 @@ double Trk::MultiComponentStateModeCalculator::findModeGlobal( const std::array<
return mode;
}
-double Trk::MultiComponentStateModeCalculator::pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double x, int i ) const
+double Trk::MultiComponentStateModeCalculator::pdf( double x, int i, std::array<std::vector< Mixture >,5>& mixture )
{
double pdf(0.);
- std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
+ auto component = mixture[i].begin();
- for ( ; component != mixtureArray[i].end() ; ++component )
+ for ( ; component != mixture[i].end() ; ++component )
pdf += component->weight * gaus( x, component->mean, component->sigma );
return pdf;
}
-double Trk::MultiComponentStateModeCalculator::d1pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double x, int i ) const
+double Trk::MultiComponentStateModeCalculator::d1pdf( double x, int i,std::array<std::vector< Mixture >,5>& mixture )
{
double result(0.);
- std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
+ auto component = mixture[i].begin();
- for( ; component != mixtureArray[i].end() ; ++component ){
+ for( ; component != mixture[i].end() ; ++component ){
double z = ( x - component->mean ) / component->sigma;
@@ -323,16 +272,14 @@ double Trk::MultiComponentStateModeCalculator::d1pdf( const std::array<std::vect
}
-double Trk::MultiComponentStateModeCalculator::d2pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
- double x, int i ) const
+double Trk::MultiComponentStateModeCalculator::d2pdf( double x, int i, std::array<std::vector< Mixture >,5>& mixture )
{
-
double result(0.);
- std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
+ auto component = mixture[i].begin();
- for( ; component != mixtureArray[i].end() ; ++component ){
+ for( ; component != mixture[i].end() ; ++component ){
double z = ( x - component->mean ) / component->sigma;
@@ -344,7 +291,7 @@ double Trk::MultiComponentStateModeCalculator::d2pdf( const std::array<std::vect
}
-double Trk::MultiComponentStateModeCalculator::gaus( double x, double mean, double sigma ) const
+double Trk::MultiComponentStateModeCalculator::gaus( double x, double mean, double sigma )
{
double normalisation = 1. / sqrt( 2. * M_PI );
@@ -359,15 +306,14 @@ double Trk::MultiComponentStateModeCalculator::gaus( double x, double mean, doub
}
-double Trk::MultiComponentStateModeCalculator::width( const std::array<std::vector< Mixture >,5>& mixtureArray ,
- int i ) const
+double Trk::MultiComponentStateModeCalculator::width( int i, std::array<std::vector< Mixture >,5>& mixture )
{
double pdf(0.);
- std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
+ auto component = mixture[i].begin();
- for ( ; component != mixtureArray[i].end() ; ++component )
+ for ( ; component != mixture[i].end() ; ++component )
pdf += component->weight * component->sigma;
return pdf;
@@ -376,24 +322,23 @@ double Trk::MultiComponentStateModeCalculator::width( const std::array<std::vect
-double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::vector< Mixture >,5>& mixtureArray,
- double &result, double xlo,
- double xhi, double value,
- double i) const
+double Trk::MultiComponentStateModeCalculator::findRoot(double &result, double xlo, double xhi, double value, double i, std::array<std::vector< Mixture >,5>& mixture, MsgStream &log)
{
// Do the root finding using the Brent-Decker method. Returns a boolean status and
// loads 'result' with our best guess at the root if true.
// Prints a warning if the initial interval does not bracket a single
// root or if the root is not found after a fixed number of iterations.
-
+
double a(xlo),b(xhi);
- double fa= pdf(mixtureArray,a,i) - value;
- double fb= pdf(mixtureArray,b,i) - value;
-
+ double fa= pdf(a,i,mixture) - value;
+ double fb= pdf(b,i,mixture) - value;
+
if(fb*fa > 0) {
- ATH_MSG_DEBUG( "BrentRootFinder::findRoot: initial interval does not bracket a root: ("
- << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb );
+
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "BrentRootFinder::findRoot: initial interval does not bracket a root: ("
+ << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endmsg;
return false;
}
@@ -402,7 +347,7 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
double c(0),d(0),e(0);
int MaxIterations = 20;
double tolerance = 1.e-6;
-
+
for(int iter= 0; iter <= MaxIterations; iter++) {
if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
@@ -413,7 +358,7 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
d = b - a;
e = b - a;
}
-
+
if (fabs (fc) < fabs (fb)) {
ac_equal = true;
a = b;
@@ -429,12 +374,14 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
if (fb == 0 || fabs(m) <= tol) {
- ATH_MSG_DEBUG( "BrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol );
- ATH_MSG_DEBUG( "BrentRootFinder: xlo = " << xlo << " < " << b << " < " << xhi );
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) ){
+ log << "BrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endmsg;
+ log << "BrentRootFinder: xlo = " << xlo << " < " << b << " < " << xhi << endmsg;
+ }
result= b;
return true;
}
-
+
if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
// Bounds decreasing too slowly: use bisection
d = m;
@@ -444,7 +391,7 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
// Attempt inverse cubic interpolation
double p, q, r;
double s = fb / fa;
-
+
if (ac_equal) {
p = 2 * m * s;
q = 1 - s;
@@ -462,7 +409,7 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
else {
p = -p;
}
-
+
double min1= 3 * m * q - fabs (tol * q);
double min2= fabs (e * q);
if (2 * p < (min1 < min2 ? min1 : min2)) {
@@ -486,11 +433,12 @@ double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::ve
else {
b += (m > 0 ? +tol : -tol);
}
- fb= pdf(mixtureArray,b, i) - value;
+ fb= pdf(b, i, mixture) - value;
}
// Return our best guess if we run out of iterations
- ATH_MSG_DEBUG( "BrentRootFinder::findRoot: maximum iterations exceeded." );
+ if( ATH_UNLIKELY( log.level() <= MSG::DEBUG ) )
+ log << "BrentRootFinder::findRoot: maximum iterations exceeded." << endmsg;
result= b;
return false;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.h b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.h
new file mode 100755
index 0000000000000000000000000000000000000000..2d73b82d51bb6df13c364f97bea171a522411368
--- /dev/null
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.h
@@ -0,0 +1,78 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***********************************************************************************
+ MultiComponentStateModeCalculator.h - description
+ ---------------------------------------------------
+begin : Thursday 6th July 2006
+author : atkinson
+email : Tom.Atkinson@cern.ch
+description : Class to calculate the mode (q/p) of a gaussian mixtureArray
+***********************************************************************************/
+
+#ifndef Trk_MultiComponentStateModeCalculator_H
+#define Trk_MultiComponentStateModeCalculator_H
+
+#include "GaudiKernel/MsgStream.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include <array>
+
+namespace Trk{
+
+class MultiComponentState;
+
+namespace MultiComponentStateModeCalculator{
+
+ struct Mixture{
+
+ // Default constructor
+ Mixture():
+ weight(0.),
+ mean(0.),
+ sigma(0.)
+ {}
+
+ // Constructor with arguments
+ Mixture( double aWeight, double aMean, double aSigma ):
+ weight( aWeight ),
+ mean( aMean ),
+ sigma( aSigma )
+ {}
+ double weight;
+ double mean;
+ double sigma;
+ };
+
+ //!< IMultiComponentStateModeCalculator interface method to calculate mode
+ Amg::VectorX calculateMode( const MultiComponentState&, MsgStream &log ) ;
+
+ //!< Private method to extract the weight, mean and sigma values from the multi-component state
+ void fillMixture( const MultiComponentState&, std::array<std::vector< Mixture >,5>& ) ;
+
+ //!< Private method to find the mode using the Newton-Raphson method based on a starting guess
+ double findMode( double, int, std::array<std::vector< Mixture >,5>& mixture, MsgStream &log ) ;
+
+ //!< Private method to determine the pdf of the cashed mixture at a given value
+ double pdf( double, int, std::array<std::vector< Mixture >,5>& mixture ) ;
+
+ //!< Private method to determine the first order derivative of the pdf at a given value
+ double d1pdf( double, int, std::array<std::vector< Mixture >,5>& mixture ) ;
+
+ //!< Private method to determine the second order derivative of the pdf at a given value
+ double d2pdf( double, int, std::array<std::vector< Mixture >,5>& mixture) ;
+
+ //!< Private method to determine the value of the a gaussian distribution at a given value
+ double gaus( double x, double mean, double sigma ) ;
+
+ double findModeGlobal(double, int,std::array<std::vector< Mixture >,5>& mixture) ;
+
+ double width( int i, std::array<std::vector< Mixture >,5>& mixture) ;
+
+ double findRoot(double &result, double xlo, double xhi, double value, double i,std::array<std::vector< Mixture >,5>& mixture, MsgStream &log) ;
+
+} // end MultiComponentStateModeCalculator namespace
+
+} // end Trk namespace
+
+#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
index 4b91eb5328b6f7e3808198fb996358344c1e995e..cb510b83fde825923be114a242e29358e4157d0b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
@@ -1,5 +1,4 @@
#include "TrkGaussianSumFilter/QuickCloseComponentsMultiStateMerger.h"
-#include "TrkGaussianSumFilter/MultiComponentStateModeCalculator.h"
#include "TrkGaussianSumFilter/MultiStateMaterialEffectsAdapter.h"
#include "TrkGaussianSumFilter/KullbackLeiblerComponentDistance.h"
#include "TrkGaussianSumFilter/CloseComponentsMultiStateMerger.h"
@@ -17,12 +16,10 @@
#include "TrkGaussianSumFilter/GsfEnergyLossUpdator.h"
#include "TrkGaussianSumFilter/GaussianSumFitter.h"
#include "TrkGaussianSumFilter/ForwardGsfFitter.h"
-#include "TrkGaussianSumFilter/GsfOutlierLogic.h"
#include "TrkGaussianSumFilter/GsfExtrapolator.h"
#include "TrkGaussianSumFilter/GsfSmoother.h"
DECLARE_COMPONENT( Trk::QuickCloseComponentsMultiStateMerger )
-DECLARE_COMPONENT( Trk::MultiComponentStateModeCalculator )
DECLARE_COMPONENT( Trk::MultiStateMaterialEffectsAdapter )
DECLARE_COMPONENT( Trk::KullbackLeiblerComponentDistance )
DECLARE_COMPONENT( Trk::CloseComponentsMultiStateMerger )
@@ -40,7 +37,5 @@ DECLARE_COMPONENT( Trk::GsfMeasurementUpdator )
DECLARE_COMPONENT( Trk::GsfEnergyLossUpdator )
DECLARE_COMPONENT( Trk::GaussianSumFitter )
DECLARE_COMPONENT( Trk::ForwardGsfFitter )
-DECLARE_COMPONENT( Trk::GsfOutlierLogic )
DECLARE_COMPONENT( Trk::GsfExtrapolator )
DECLARE_COMPONENT( Trk::GsfSmoother )
-