Skip to content
Snippets Groups Projects
Commit b8e0c1c4 authored by Adam Edward Barton's avatar Adam Edward Barton
Browse files

Merge branch 'MT_MCSModeCalculator' into 'master' 

Threadsafe  MultiComponentStateCombiner: ATLASRECTS-4612

See merge request atlas/athena!14669
parents 167329a1 d682640a
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 430 additions and 643 deletions
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
+ 16
− 17
View file @ b8e0c1c4
......@@ -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
......
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfOutlierLogic.h deleted 100755 → 0
+ 0
− 52
View file @ 167329a1
/*
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
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IGsfOutlierLogic.h deleted 100755 → 0
+ 0
− 43
View file @ 167329a1
/*
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
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateCombiner.h
+ 4
− 8
View file @ b8e0c1c4
......@@ -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
......
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h deleted 100755 → 0
+ 0
− 46
View file @ 167329a1
/*
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
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateCombiner.h
+ 9
− 18
View file @ b8e0c1c4
......@@ -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
......
Tracking/TrkFitter/TrkGaussianSumFilter/src/CloseComponentsMultiStateMerger.cxx
+ 7
− 9
View file @ b8e0c1c4
......@@ -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;
......
Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfOutlierLogic.cxx deleted 100755 → 0
+ 0
− 76
View file @ 167329a1
/*
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;
}
Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateCombiner.cxx
+ 40
− 57
View file @ b8e0c1c4
......@@ -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;
}
Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateModeCalculator.h → Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.h
+ 78
− 0
View file @ b8e0c1c4
......@@ -14,38 +14,16 @@ description : Class to calculate the mode (q/p) of a gaussian mixtureAr
#ifndef Trk_MultiComponentStateModeCalculator_H
#define Trk_MultiComponentStateModeCalculator_H
#include "AthenaBaseComps/AthAlgTool.h"
#include "TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h"
#include <atomic>
#include "GaudiKernel/MsgStream.h"
#include "EventPrimitives/EventPrimitives.h"
#include <array>
#include <vector>
namespace Trk{
class MultiComponentState;
class MultiComponentStateModeCalculator : public AthAlgTool, virtual public IMultiComponentStateModeCalculator{
namespace MultiComponentStateModeCalculator{
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
......@@ -66,53 +44,34 @@ class MultiComponentStateModeCalculator : public AthAlgTool, virtual public IMul
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( std::array<std::vector< Mixture >,5>& mixtureArray,
const MultiComponentState& ) const;
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( const std::array<std::vector< Mixture >,5>& mixtureArray,
double, int ) const;
double findMode( double, int, std::array<std::vector< Mixture >,5>& mixture, MsgStream &log ) ;
//!< 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 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( const std::array<std::vector< Mixture >,5>& mixtureArray,
double, int ) const;
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( const std::array<std::vector< Mixture >,5>& mixtureArray,
double, int) const;
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 ) 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;
};
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
......
Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
+ 0
− 5
View file @ b8e0c1c4
#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 )
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment