diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/ElectronEoverPTracking.py b/InnerDetector/InDetExample/InDetAlignExample/share/ElectronEoverPTracking.py
index 83083ddcf3278618a6326c2a9cd1d030a71e2996..e8fcde204bfe9782014ef50627904f53493c8157 100755
--- a/InnerDetector/InDetExample/InDetAlignExample/share/ElectronEoverPTracking.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/ElectronEoverPTracking.py
@@ -168,14 +168,6 @@ GsfExtrapolator = Trk__GsfExtrapolator(name = 'GsfExtra
ToolSvc += GsfExtrapolator
-from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__BremFind
-BremFind = Trk__BremFind(name = 'BremFind',
- UseCalibration = True,
- ValidationMode = True )
-
-
-ToolSvc += BremFind
-
from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__GaussianSumFitter
GSFTrackFitter = Trk__GaussianSumFitter(name = 'GSFTrackFitter',
ToolForExtrapolation = GsfExtrapolator,
@@ -186,8 +178,6 @@ GSFTrackFitter = Trk__GaussianSumFitter(name = 'GSFTrackFitte
RefitOnMeasurementBase = True,
DoHitSorting = True,
ValidationMode = False,
- BremFind = BremFind,
- runBremFinder = False,
OutputLevel = 3)
# --- end of fitter loading
ToolSvc += GSFTrackFitter
diff --git a/InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/share/ElectronEoverPTracking.py b/InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/share/ElectronEoverPTracking.py
index 9de3dd15492b83c611bf3e21069c5764f60dab59..50868637b9b5649fc379fbf7683f801fd721cbea 100755
--- a/InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/share/ElectronEoverPTracking.py
+++ b/InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/share/ElectronEoverPTracking.py
@@ -170,14 +170,6 @@ GsfExtrapolator = Trk__GsfExtrapolator(name = 'GsfExtra
ToolSvc += GsfExtrapolator
-from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__BremFind
-BremFind = Trk__BremFind(name = 'BremFind',
- UseCalibration = True,
- ValidationMode = True )
-
-
-ToolSvc += BremFind
-
from TrkGaussianSumFilter.TrkGaussianSumFilterConf import Trk__GaussianSumFitter
GSFTrackFitter = Trk__GaussianSumFitter(name = 'GSFTrackFitter',
ToolForExtrapolation = GsfExtrapolator,
@@ -188,8 +180,6 @@ GSFTrackFitter = Trk__GaussianSumFitter(name = 'GSFTrackFitte
RefitOnMeasurementBase = True,
DoHitSorting = True,
ValidationMode = False,
- BremFind = BremFind,
- runBremFinder = False,
OutputLevel = 3)
# --- end of fitter loading
ToolSvc += GSFTrackFitter
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/BremFind.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/BremFind.h
deleted file mode 100644
index 5d549bab03fb4a9454c573798dded0ca9bdefe04..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/BremFind.h
+++ /dev/null
@@ -1,229 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/********************************************************************************
- BremFind.h * description
-
-completed : Tuesday 18th November 2008
-author : Tony Shao
-email : Qi.Tao.Shao@cern.ch
-description : Class for finding brem points in the inner detector using the GSF
-**********************************************************************************/
-
-#ifndef Trk_BremFind_H
-#define Trk_BremFind_H
-
-#include "TrkGaussianSumFilter/IBremsstrahlungFinder.h"
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "GaudiKernel/ServiceHandle.h"
-
-#include "AthContainers/DataVector.h"
-#include "TrkFitterUtils/FitterTypes.h"
-#include "TrkParameters/TrackParameters.h"
-
-#include "TrkExInterfaces/IPropagator.h"
-#include "TrkDetDescrInterfaces/ITrackingGeometrySvc.h"
-#include "TrkGeometry/MagneticFieldProperties.h"
-
-#include "GeoPrimitives/GeoPrimitives.h"
-
-#include "StoreGate/ReadHandleKey.h"
-#include "xAODEventInfo/EventInfo.h"
-
-#include <vector>
-
-//class ITrackingGeometrySvc;
-class TTree;
-
-namespace Trk {
-
-class IMultiComponentStateCombiner;
-class MultiComponentStateOnSurface;
-class TrackStateOnSurface;
-class EstimatedBremOnTrack;
-class Layer;
-class TrackingVolume;
-class TrackingGeometry;
-class EntryLayerProvider;
-class AbstractVolume;
-class Surface;
-class QoverPBremFit;
-
-struct GraphParameter {
- double constant;
- std::vector<double> coefficient{};
- std::vector<double> value{};
- std::vector<double> width{};
-};
-
-struct Element {
- double value{};
- double width{};
- bool sign{};
-};
-
-
-#define TRKFGSF_VALSURFACES 100
-#define TRKGSF_VALSTATES 24
-#define POTENTIAL_BREM_THRESHOLD 0.0001
-#define MIN_KINK_THRESHOLD 0.005
-#define BREM_AMPLITUDE 1.0 //in GeV
-#define BREM_SEPARATION 0.005
-#define COEFFICIENT_THRESHOLD 0.005
-#define BOUND_THRESHOLD 20.0 //in percentage
-
-class BremFind : public AthAlgTool, virtual public IBremsstrahlungFinder{
- public:
-
- BremFind(const std::string&, const std::string&, const IInterface*);
-
- virtual ~BremFind () {};
-
- StatusCode initialize();
-
- /** AlgTool finalise method */
- StatusCode finalize();
-
- //Find the brem points from a forward and smooth track fit
- virtual void BremFinder(const Trk::ForwardTrajectory&, const Trk::SmoothedTrajectory&);
-
- virtual void BremFinder(const Trk::ForwardTrajectory&, const Trk::SmoothedTrajectory&, bool);
-
- virtual void CombineParameters();
-
- virtual void GetBremLocations();
-
- virtual double GetBound(double, int);
-
- virtual double BremSignificance(double, QoverPBremFit*, QoverPBremFit*);
-
- virtual double GetPhiFromValue(double, QoverPBremFit*);
-
- virtual int ClosestMeasurement(double, QoverPBremFit*);
-
- virtual void CalibrateValue();
-
- virtual const Surface* ClosestSurface(double, double, double, const TrackParameters&);
-
- virtual Amg::Vector3D SurfacePosition(const TrackParameters&, const Surface&, double, double, double);
-
- virtual StatusCode retrieveTrackingGeometry();
-
- virtual int GetNBrems();
-
- virtual void MakeTrackStateOnSurfaces ();
-
- virtual const TrackStateOnSurface* GetEstimatedBremOnTrack(int);
-
-
- private:
-
- double Graph_value(GraphParameter, double);
-
- double GetMaximumX(GraphParameter, double, double);
-
- double GetMaximumY(GraphParameter, double, double);
-
- double UniqueAngle(double);
-
- void SeparationQuality(std::vector<double>*, std::vector<double>*);
-
- std::pair<double,double> ProbabilityScore(double, double, double, double);
-
- //defining the tools used in this package
- ToolHandle<IMultiComponentStateCombiner> m_stateCombiner;
-
- ToolHandle<IPropagator> m_propagator;
- MagneticFieldProperties m_fieldProperties;
-
- ServiceHandle<ITrackingGeometrySvc> m_trackingGeometrySvc; //!< Name of the TrackingGeometrySvc
- std::string m_trackingGeometryName; //!< TrackingGeometry to be retrieved from DetectorStore
- mutable const TrackingGeometry* m_trackingGeometry; //!< The underlying TrackingGeometry
-
- //Validation variables
-
- bool m_validationMode;
- std::string m_validationTreeName;
- std::string m_validationTreeName2;
- std::string m_validationTreeDescription;
- std::string m_validationTreeFolder;
- std::string m_validationTreeFolder2;
-
- TTree* m_validationTree;
- TTree* m_validationTree2;
-
- //Calibration option
- bool m_useCalibration;
-
- //Propagation option
- bool m_usePropagate;
-
- //Class variables
- GraphParameter m_forwardparameters;
- GraphParameter m_smoothedparameters;
- GraphParameter m_combinedparameters;
-
- double m_perigee_1overP;
- double m_perigee_Phi;
- double m_perigee_Theta;
- double m_perigee_d0;
- double m_perigee_z0;
-
- std::vector<double> *m_brem_value; //could be R or Z;
- std::vector<double> *m_brem_phi;
- std::vector<double> *m_brem_theta;
- std::vector<double> *m_brem_energy;
- std::vector<double> *m_brem_UpperBound;
- std::vector<double> *m_brem_LowerBound;
- std::vector<double> *m_forward_kink;
- std::vector<double> *m_smoothed_kink;
- std::vector<double> *m_brem_significance;
- std::vector<double> *m_brem_valueCalibrated;
- std::vector<double> *m_surfaceX;
- std::vector<double> *m_surfaceY;
- std::vector<double> *m_surfaceZ;
-
-
- int m_nBrems;
- bool m_Z_mode;
-
- std::vector<const Trk::TrackStateOnSurface*> m_brem_trackStateOnSurface;
- std::vector<const Trk::TrackParameters*> m_brem_TrackParameters;
- std::vector<const Trk::EstimatedBremOnTrack*> m_EstimatedBremOnTrack;
- std::vector<const Trk::Surface*> m_brem_surfaces;
- std::vector<const Trk::Layer*> m_brem_layers;
- std::vector<Amg::Vector3D> m_surface_positions;
-
- int m_event_ID;
- double m_forwardparameter_constant;
- std::vector<double> *m_forwardparameter_coefficient;
- std::vector<double> *m_forwardparameter_value;
- std::vector<double> *m_forwardparameter_width;
- double m_smoothparameter_constant;
- std::vector<double> *m_smoothparameter_coefficient;
- std::vector<double> *m_smoothparameter_value;
- std::vector<double> *m_smoothparameter_width;
-
- std::vector<double> *m_forward_1overP;
- std::vector<double> *m_forward_1overPerr;
- std::vector<double> *m_forward_value;
- std::vector<double> *m_smooth_1overP;
- std::vector<double> *m_smooth_1overPerr;
- std::vector<double> *m_smooth_value;
-
- std::vector<double> *m_KinkSeparationScores;
- std::vector<double> *m_KinkSeparationScoresErr;
-
- QoverPBremFit *m_forwardBremFit;
- QoverPBremFit *m_smoothedBremFit;
-
- SG::ReadHandleKey<xAOD::EventInfo> m_readKey;
-
- };
-
-}
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
index 4665c9212f6533f603d480b5bcfbb89c742f222c..9c80d601ac973a40df7b2e04264c919dd0c364d6 100755
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
@@ -40,7 +40,6 @@ class TrackFitInputPreparator;
class IForwardGsfFitter;
class IGsfSmoother;
class Track;
-class IBremsstrahlungFinder;
#define TRKFGSF_VALSURFACES 100
#define TRKGSF_VALSTATES 24
@@ -139,19 +138,12 @@ class GaussianSumFitter : virtual public ITrackFitter, public AthAlgTool {
bool m_makePerigee;
PropDirection m_directionToPerigee;
-
- bool m_runBremFinder;
-
bool m_refitOnMeasurementBase;
-
bool m_doHitSorting;
TrkParametersComparisonFunction* m_trkParametersComparisonFunction;
std::vector<double> m_sortingReferencePoint;
ToolHandle<IMultiComponentStateCombiner> m_stateCombiner;
- ToolHandle<IBremsstrahlungFinder> m_BremFind;
- ToolHandle<IBremsstrahlungFinder> m_BremFind2;
-
ServiceHandle<IChronoStatSvc> m_chronoSvc;
TrackFitInputPreparator* m_inputPreparator;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBremsstrahlungFinder.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBremsstrahlungFinder.h
deleted file mode 100644
index c54728fcfd8412c849cd66b5bafe70752ff535d9..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBremsstrahlungFinder.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/*********************************************************************************
- IBremsstrahlungFinder.h - description
- ---------------------------------------------
-created : Thursday 8th January 2009
-author : amorley
-email : Anthony.Morley@cern.ch
-description : Interface for bremsstrahlung finding using the GSF
-*********************************************************************************/
-
-
-#ifndef Trk_IBremsstrahlungFinder_H
-#define Trk_IBremsstrahlungFinder_H
-
-#include "GaudiKernel/IAlgTool.h"
-
-#include "TrkParameters/TrackParameters.h"
-#include "TrkFitterUtils/FitterTypes.h"
-
-#include "GeoPrimitives/GeoPrimitives.h"
-
-namespace Trk{
-
-class MultiComponentState;
-class TrackStateOnSurface;
-class Surface;
-class QoverPBremFit;
-
-
-//!< Interface ID for IBremsstrahlungFinder
-static const InterfaceID IID_IBremsstrahlungFinder( "IBremsstrahlungFinder", 1, 0 );
-
-class IBremsstrahlungFinder : virtual public IAlgTool{
-
- public:
-
- //!< Virtual destructor
- virtual ~IBremsstrahlungFinder(){};
-
- //!< IAlgTool-AlgTool interface
- static const InterfaceID& interfaceID() { return IID_IBremsstrahlungFinder; };
-
- //Find the brem points from a forward and smooth track fit
- virtual void BremFinder(const Trk::ForwardTrajectory&, const Trk::SmoothedTrajectory&) = 0;
-
- virtual void BremFinder(const Trk::ForwardTrajectory&, const Trk::SmoothedTrajectory&, bool) = 0;
-
- virtual void CombineParameters() = 0;
-
- virtual void GetBremLocations() = 0;
-
- virtual double GetBound(double, int) = 0;
-
- virtual double BremSignificance(double, QoverPBremFit*, QoverPBremFit*) = 0;
-
- virtual double GetPhiFromValue(double, QoverPBremFit*) = 0;
-
- virtual int ClosestMeasurement(double, QoverPBremFit*) = 0;
-
- virtual void CalibrateValue() = 0;
-
- virtual const Surface* ClosestSurface(double, double, double, const TrackParameters&) = 0;
-
- virtual Amg::Vector3D SurfacePosition(const TrackParameters&, const Surface&, double, double, double) = 0;
-
- virtual StatusCode retrieveTrackingGeometry() = 0;
-
- virtual int GetNBrems() = 0;
-
- virtual void MakeTrackStateOnSurfaces () = 0;
-
- virtual const TrackStateOnSurface* GetEstimatedBremOnTrack(int) = 0;
-
-};
-
-} // end Trk namespace
-
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/BremFind.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/BremFind.cxx
deleted file mode 100644
index 88733889ac8814888cee8ebf807088822d53bcb0..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/BremFind.cxx
+++ /dev/null
@@ -1,1001 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/********************************************************************************
- BremFind.cxx * description
-
-completed : Tuesday 18th November 2008
-author : Tony Shao
-email : Qi.Tao.Shao@cern.ch
-description : Class for finding brem points in the inner detector using the GSF
-**********************************************************************************/
-
-
-
-#include "TrkGaussianSumFilter/BremFind.h"
-#include "QoverPBremFit.h"
-
-#include <math.h>
-#include "TMath.h"
-
-#include "TrkGaussianSumFilter/IMultiComponentStateCombiner.h"
-
-#include "TrkMultiComponentStateOnSurface/MultiComponentStateOnSurface.h"
-#include "TrkParameters/TrackParameters.h"
-
-#include "TrkTrack/Track.h"
-#include "TrkMaterialOnTrack/EstimatedBremOnTrack.h"
-
-#include "TrkGeometry/Layer.h"
-#include "TrkGeometry/TrackingVolume.h"
-#include "TrkGeometry/TrackingGeometry.h"
-#include "TrkGeometry/MagneticFieldProperties.h"
-
-#include "TrkVolumes/AbstractVolume.h"
-#include "TrkSurfaces/Surface.h"
-#include "TrkSurfaces/CylinderSurface.h"
-
-#include "GaudiKernel/ITHistSvc.h"
-
-#include "TTree.h"
-
-
-
-Trk::BremFind::BremFind(const std::string& type, const std::string& name, const IInterface* parent):
- AthAlgTool(type, name, parent),
- m_stateCombiner("Trk::MultiComponentStateCombiner"),
- m_propagator("Trk::IntersectorWrapper/IntersectorWrapper"),
- //m_fieldProperties
- m_trackingGeometrySvc("TrackingGeometrySvc","AtlasTrackingGeometrySvc"),
- m_trackingGeometryName("AtlasTrackingGeometry"),
- m_trackingGeometry(nullptr),
- //m_validationMode
- m_validationTreeName("BremInfo"),
- m_validationTreeName2("BremInfoZ"),
- m_validationTreeDescription("Brem Information"),
- m_validationTreeFolder("/valGSF2/BremInfo"),
- m_validationTreeFolder2("/valGSF2/BremInfoZ"),
- m_validationTree(nullptr),
- m_validationTree2(nullptr),
- //m_useCalibration{},
- //m_usePropagate{},
- m_forwardparameters{},
- m_smoothedparameters{},
- m_combinedparameters{},
- m_perigee_1overP{},
- m_perigee_Phi{},
- m_perigee_Theta{},
- m_perigee_d0{},
- m_perigee_z0{},
- m_brem_value{},
- m_brem_phi{},
- m_brem_theta{},
- m_brem_energy{},
- m_brem_UpperBound{},
- m_brem_LowerBound{},
- m_forward_kink{},
- m_smoothed_kink{},
- m_brem_significance{},
- m_brem_valueCalibrated{},
- m_surfaceX{},
- m_surfaceY{},
- m_surfaceZ{},
- m_nBrems{},
- m_Z_mode{},
- //vectors could be initialised here
- m_event_ID{},
- m_forwardparameter_constant{},
- m_forwardparameter_coefficient{},
- m_forwardparameter_value{},
- m_forwardparameter_width{},
- m_smoothparameter_constant{},
- m_smoothparameter_coefficient{},
- m_smoothparameter_value{},
- m_smoothparameter_width{},
- m_forward_1overP{},
- m_forward_1overPerr{},
- m_forward_value{},
- m_smooth_1overP{},
- m_smooth_1overPerr{},
- m_smooth_value{},
- m_KinkSeparationScores{},
- m_KinkSeparationScoresErr{},
- m_forwardBremFit{},
- m_smoothedBremFit{}
-{
- declareInterface<IBremsstrahlungFinder>(this);
- //jobOptions Variables
- declareProperty("StateCombiner", m_stateCombiner );
- declareProperty("TrackingGeometrySvc", m_trackingGeometrySvc);
- declareProperty("Propagator", m_propagator);
- declareProperty("UseCalibration", m_useCalibration=true);
- declareProperty("ValidationMode", m_validationMode=false);
- declareProperty("UseSurfacePropagation", m_usePropagate=false);
- declareProperty("TreeFolderName", m_validationTreeFolder);
- declareProperty("TreeFolderName2", m_validationTreeFolder2);
- declareProperty("EventInfoKey", m_readKey = "EventInfo");
-}
-
-
-StatusCode Trk::BremFind::initialize()
-{
- // The TrackingGeometrySvc ------------------------------------------------------
- ATH_CHECK (m_trackingGeometrySvc.retrieve());
-
- msg(MSG::DEBUG) << "Retrieved service " << m_trackingGeometrySvc << endmsg;
- m_trackingGeometryName = m_trackingGeometrySvc->trackingGeometryName();
-
-
- // Request the state combiner
- ATH_CHECK ( m_stateCombiner.retrieve() );
-
- //Retrieve the propagator
- ATH_CHECK (m_propagator.retrieve());
-
- m_fieldProperties= Trk::MagneticFieldProperties(Trk::FullField);
-
-
-
- //The Validation step ----------------------------------------------------------
-
- if (m_validationMode) {
-
- if (not m_validationTree) {
- //Crate a new tree if there doesn't exist one already
- m_validationTree = new TTree( m_validationTreeName.c_str(), m_validationTreeDescription.c_str() );
-
-
- m_validationTree->Branch("nBrems" , &m_nBrems, "nBrems/I");
- m_validationTree->Branch("Brem_Radii", &m_brem_value);
- m_validationTree->Branch("Brem_Phi", &m_brem_phi);
- m_validationTree->Branch("Brem_Theta", &m_brem_theta);
- m_validationTree->Branch("Brem_Energy", &m_brem_energy);
- m_validationTree->Branch("Brem_UpperBound", &m_brem_UpperBound);
- m_validationTree->Branch("Brem_LowerBound", &m_brem_LowerBound);
- m_validationTree->Branch("Forward_Kink", &m_forward_kink);
- m_validationTree->Branch("Backward_Kink", &m_smoothed_kink);
- m_validationTree->Branch("Brem_Significance", &m_brem_significance);
- m_validationTree->Branch("Brem_RadiiCalibrated", &m_brem_valueCalibrated);
- m_validationTree->Branch("SurfaceX", &m_surfaceX);
- m_validationTree->Branch("SurfaceY", &m_surfaceY);
- m_validationTree->Branch("SurfaceZ", &m_surfaceZ);
- m_validationTree->Branch("EventID" , &m_event_ID, "EventID/I");
- m_validationTree->Branch("ForwardParameterConstant", &m_forwardparameter_constant,"ForwardParameterConstant/D");
- m_validationTree->Branch("SmoothParameterConstant", &m_smoothparameter_constant,"SmoothParameterConstant/D");
- m_validationTree->Branch("ForwardParameterCoefficient", &m_forwardparameter_coefficient);
- m_validationTree->Branch("ForwardParameterRadii", &m_forwardparameter_value);
- m_validationTree->Branch("ForwardParameterWidth", &m_forwardparameter_width);
- m_validationTree->Branch("SmoothParameterCoefficient", &m_smoothparameter_coefficient);
- m_validationTree->Branch("SmoothParameterRadii", &m_smoothparameter_value);
- m_validationTree->Branch("SmoothParameterWidth", &m_smoothparameter_width);
- m_validationTree->Branch("Forward1overP", &m_forward_1overP);
- m_validationTree->Branch("Forward1overPerr", &m_forward_1overPerr);
- m_validationTree->Branch("ForwardRadii", &m_forward_value);
- m_validationTree->Branch("Smooth1overP", &m_smooth_1overP);
- m_validationTree->Branch("Smooth1overPerr", &m_smooth_1overPerr);
- m_validationTree->Branch("SmoothRadii", &m_smooth_value);
- m_validationTree->Branch("SeparationQuality", &m_KinkSeparationScores);
- m_validationTree->Branch("SeparationQualityErr", &m_KinkSeparationScoresErr);
-
- }
-
-
-
-
- if (not m_validationTree2) {
- //Crate a new tree if there doesn't exist one already
- m_validationTree2 = new TTree( m_validationTreeName2.c_str(), m_validationTreeDescription.c_str() );
-
-
- m_validationTree2->Branch("nBrems" , &m_nBrems, "nBrems/I");
- m_validationTree2->Branch("Brem_Z", &m_brem_value);
- m_validationTree2->Branch("Brem_Phi", &m_brem_phi);
- m_validationTree2->Branch("Brem_Theta", &m_brem_theta);
- m_validationTree2->Branch("Brem_Energy", &m_brem_energy);
- m_validationTree2->Branch("Brem_UpperBound", &m_brem_UpperBound);
- m_validationTree2->Branch("Brem_LowerBound", &m_brem_LowerBound);
- m_validationTree2->Branch("Forward_Kink", &m_forward_kink);
- m_validationTree2->Branch("Backward_Kink", &m_smoothed_kink);
- m_validationTree2->Branch("Brem_Significance", &m_brem_significance);
- m_validationTree2->Branch("Brem_ZCalibrated", &m_brem_valueCalibrated);
- m_validationTree2->Branch("SurfaceX", &m_surfaceX);
- m_validationTree2->Branch("SurfaceY", &m_surfaceY);
- m_validationTree2->Branch("SurfaceZ", &m_surfaceZ);
- m_validationTree2->Branch("EventID" , &m_event_ID, "EventID/I");
- m_validationTree2->Branch("ForwardParameterConstant", &m_forwardparameter_constant,"ForwardParameterConstant/D");
- m_validationTree2->Branch("SmoothParameterConstant", &m_smoothparameter_constant,"SmoothParameterConstant/D");
- m_validationTree2->Branch("ForwardParameterCoefficient", &m_forwardparameter_coefficient);
- m_validationTree2->Branch("ForwardParameterZ", &m_forwardparameter_value);
- m_validationTree2->Branch("ForwardParameterWidth", &m_forwardparameter_width);
- m_validationTree2->Branch("SmoothParameterCoefficient", &m_smoothparameter_coefficient);
- m_validationTree2->Branch("SmoothParameterZ", &m_smoothparameter_value);
- m_validationTree2->Branch("SmoothParameterWidth", &m_smoothparameter_width);
- m_validationTree2->Branch("Forward1overP", &m_forward_1overP);
- m_validationTree2->Branch("Forward1overPerr", &m_forward_1overPerr);
- m_validationTree2->Branch("ForwardZ", &m_forward_value);
- m_validationTree2->Branch("Smooth1overP", &m_smooth_1overP);
- m_validationTree2->Branch("Smooth1overPerr", &m_smooth_1overPerr);
- m_validationTree2->Branch("SmoothZ", &m_smooth_value);
- m_validationTree2->Branch("SeparationQuality", &m_KinkSeparationScores);
- m_validationTree2->Branch("SeparationQualityErr", &m_KinkSeparationScoresErr);
-
- }
-
-
- //register the Tree
- ITHistSvc* tHistSvc = 0;
- if (service("THistSvc",tHistSvc).isFailure()) {
- msg(MSG::ERROR) << "initialize() Could not find Hist Service -> Switching ValidationMode Off !" << endmsg;
- delete m_validationTree;
- delete m_validationTree2;
- m_validationTree=0;
- m_validationTree2=0;
- }
- if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
- msg(MSG::ERROR)<<"initialize() Could not register the validation Tree -> Switching ValidationMode off !" << endmsg;
- delete m_validationTree;
- m_validationTree = 0;
- }
- if ((tHistSvc->regTree(m_validationTreeFolder2, m_validationTree2)).isFailure()) {
- msg(MSG::ERROR)<<"initialize() Could not register the validation Tree 2 -> Switching ValidationMode off !" << endmsg;
- delete m_validationTree2;
- m_validationTree2 = 0;
- }
-
- }
-
- //---------------------------- end of validation mode ------------------------------------
-
-
- ATH_CHECK( m_readKey.initialize() );
-
- ATH_MSG_DEBUG( "Initialisation of " << name() << " was successful" );
-
- return StatusCode::SUCCESS;
-
-}
-
-
-StatusCode Trk::BremFind::finalize(){
-
- ATH_MSG_DEBUG( "Finalisation of " << name() << " was successful" );
-
- return StatusCode::SUCCESS;
-
-}
-
-void Trk::BremFind::BremFinder(const Trk::ForwardTrajectory& forwardTrajectory, const Trk::SmoothedTrajectory& smoothedTrajectory) {
- BremFinder(forwardTrajectory, smoothedTrajectory, false);
-}
-
-void Trk::BremFind::BremFinder(const Trk::ForwardTrajectory& forwardTrajectory, const Trk::SmoothedTrajectory& smoothedTrajectory, bool Z_mode)
-{
- m_Z_mode = Z_mode;
-
- if (!m_trackingGeometry)
- retrieveTrackingGeometry();
-
-
- msg(MSG::DEBUG) << "Entering BremFinder Function " << endmsg;
-
- m_brem_value = new std::vector<double>;
- m_brem_phi = new std::vector<double>;
- m_brem_theta = new std::vector<double>;
- m_brem_energy = new std::vector<double>;
- m_brem_UpperBound = new std::vector<double>;
- m_brem_LowerBound = new std::vector<double>;
- m_forward_kink = new std::vector<double>;
- m_smoothed_kink = new std::vector<double>;
- m_brem_significance = new std::vector<double>;
- m_brem_valueCalibrated = new std::vector<double>;
- m_surfaceX = new std::vector<double>;
- m_surfaceY = new std::vector<double>;
- m_surfaceZ = new std::vector<double>;
- m_brem_trackStateOnSurface.clear();
- m_EstimatedBremOnTrack.clear();
- m_brem_TrackParameters.clear();
- m_brem_surfaces.clear();
- m_brem_layers.clear();
- m_surface_positions.clear();
-
- m_forwardparameter_coefficient = new std::vector<double>;
- m_forwardparameter_value = new std::vector<double>;
- m_forwardparameter_width = new std::vector<double>;
- m_smoothparameter_coefficient = new std::vector<double>;
- m_smoothparameter_value = new std::vector<double>;
- m_smoothparameter_width = new std::vector<double>;
- m_forward_1overP = new std::vector<double>;
- m_forward_1overPerr = new std::vector<double>;
- m_forward_value = new std::vector<double>;
- m_smooth_1overP = new std::vector<double>;
- m_smooth_1overPerr = new std::vector<double>;
- m_smooth_value = new std::vector<double>;
- m_KinkSeparationScores = new std::vector<double>;
- m_KinkSeparationScoresErr = new std::vector<double>;
-
-
- //Obtain the parameters of the Tanh fit for both forward and smooth trajectories
- m_forwardBremFit = new QoverPBremFit(1,m_Z_mode);
- m_forwardparameters = m_forwardBremFit->GetParameters(forwardTrajectory);
- m_smoothedBremFit = new QoverPBremFit(-1,m_Z_mode);
- m_smoothedparameters = m_smoothedBremFit->GetParameters(smoothedTrajectory);
- //GetPerigee functions will only make sense when using the smoothed trajectory
- m_perigee_1overP = m_smoothedBremFit->GetPerigee_1overP();
- m_perigee_Phi = m_smoothedBremFit->GetPerigee_Phi();
- m_perigee_Theta = m_smoothedBremFit->GetPerigee_Theta();
- m_perigee_d0 = m_smoothedBremFit->GetPerigee_d0();
- m_perigee_z0 = m_smoothedBremFit->GetPerigee_z0();
-
-
- CombineParameters();
-
- GetBremLocations();
-
- m_nBrems = m_brem_value->size();
-
- for (int bound_counter(0); bound_counter < (int) m_brem_value->size(); bound_counter++) {
- m_brem_UpperBound->push_back(GetBound((*m_brem_value)[bound_counter],1));
- m_brem_LowerBound->push_back(GetBound((*m_brem_value)[bound_counter],-1));
- m_brem_significance->push_back(BremSignificance((*m_brem_value)[bound_counter],m_forwardBremFit,m_smoothedBremFit));
- }
-
- //Calibration of the brem value
- if (m_useCalibration) {
- CalibrateValue();
- }
- else {
- for (int copy_counter(0); copy_counter < (int) m_brem_value->size(); copy_counter++) {
- m_brem_valueCalibrated->push_back((*m_brem_value)[copy_counter]);
- }
- }
-
- for (int fill_counter(0); fill_counter < (int) m_brem_value->size(); fill_counter++) {
- m_brem_theta->push_back(m_perigee_Theta);
- m_brem_phi->push_back(GetPhiFromValue((*m_brem_valueCalibrated)[fill_counter],m_smoothedBremFit));
- }
-
- //double closest_floorF,closest_floorS;
- //int closestSurfaceNoF(0);
- int closestSurfaceNoS(0);
-
- for (int brem_counter2(0); brem_counter2 < (int) m_brem_value->size(); brem_counter2++) {
-
- //The closest measurement numbers to the brem location
- //closestSurfaceNoF=ClosestMeasurement((*m_brem_valueCalibrated)[brem_counter2],m_forwardBremFit);
- closestSurfaceNoS=ClosestMeasurement((*m_brem_valueCalibrated)[brem_counter2],m_smoothedBremFit);
-
- //closest_floorF = m_forwardBremFit->GetValue(closestSurfaceNoF);
- //closest_floorS = m_smoothedBremFit->GetValue(closestSurfaceNoS);
-
- //Filling the associated surfaces for each brem
- const Trk::TrackParameters* closestTrackParameters = m_smoothedBremFit->GetTrackParameters(closestSurfaceNoS);
- m_brem_surfaces.push_back(ClosestSurface((*m_brem_value)[brem_counter2],(*m_brem_phi)[brem_counter2],(*m_brem_theta)[brem_counter2],*closestTrackParameters));
-
- m_surface_positions.push_back(SurfacePosition(*closestTrackParameters,*m_brem_surfaces[brem_counter2],(*m_brem_value)[brem_counter2],(*m_brem_phi)[brem_counter2],(*m_brem_theta)[brem_counter2]));
-
-
- m_surfaceX->push_back(m_surface_positions[brem_counter2].x());
- m_surfaceY->push_back(m_surface_positions[brem_counter2].y());
- m_surfaceZ->push_back(m_surface_positions[brem_counter2].z());
-
-
- }
-
- double remaining_energy(1/m_perigee_1overP);
-
-
- //This loop fills the EstimatedBremOnTrack vector
-
- for (int brem_counter(0); brem_counter < (int) m_brem_value->size(); brem_counter++) {
-
-
- Trk::EstimatedBremOnTrack *estimatedBremOnTrack(0);
- remaining_energy -= (*m_brem_energy)[brem_counter];
- double retainedEnFraction(remaining_energy*m_perigee_1overP);
-
- double sigma1overPsquared(m_smoothedBremFit->Get1overPerror(0));
-
-
- if (m_brem_layers[brem_counter] && m_brem_surfaces[brem_counter] && m_brem_TrackParameters[brem_counter]) {
- estimatedBremOnTrack = new EstimatedBremOnTrack(m_brem_layers[brem_counter]->thickness(),
- retainedEnFraction,0.0,
- sigma1overPsquared,
- //*m_brem_surfaces[brem_counter],
- *m_brem_TrackParameters[brem_counter]->associatedSurface().clone(),
- bothway);
-
- }
-
- m_EstimatedBremOnTrack.push_back(estimatedBremOnTrack);
-
- }
-
- //Creates a vector of TrackStateOnSurface
- MakeTrackStateOnSurfaces();
-
-
- //* Retrieve the event info for later syncrinization
- SG::ReadHandle< xAOD::EventInfo> eventInfo (m_readKey);
- if (!eventInfo.isValid()) {
- msg(MSG::ERROR) << "Could not retrieve event info" << endmsg;
- }
-
- m_event_ID = eventInfo->eventNumber();
-
- //Fill the TanH coefficients and graph values
- for (int forward_fill(0); forward_fill < (int) m_forwardparameters.coefficient.size(); forward_fill++) {
- m_forwardparameter_coefficient->push_back(m_forwardparameters.coefficient[forward_fill]);
- m_forwardparameter_value->push_back(m_forwardparameters.value[forward_fill]);
- m_forwardparameter_width->push_back(m_forwardparameters.width[forward_fill]);
- }
-
- for (int smooth_fill(0); smooth_fill < (int) m_smoothedparameters.coefficient.size(); smooth_fill++) {
- m_smoothparameter_coefficient->push_back(m_smoothedparameters.coefficient[smooth_fill]);
- m_smoothparameter_value->push_back(m_smoothedparameters.value[smooth_fill]);
- m_smoothparameter_width->push_back(m_smoothedparameters.width[smooth_fill]);
- }
-
- for (int forwardgraph_fill(0); forwardgraph_fill < m_forwardBremFit->GetnSurfaces(); forwardgraph_fill++) {
- m_forward_1overP->push_back(m_forwardBremFit->Get1overP(forwardgraph_fill));
- m_forward_1overPerr->push_back(m_forwardBremFit->Get1overPerror(forwardgraph_fill));
- m_forward_value->push_back(m_forwardBremFit->GetValue(forwardgraph_fill));
- }
-
- for (int smoothgraph_fill(0); smoothgraph_fill < m_smoothedBremFit->GetnSurfaces(); smoothgraph_fill++) {
- m_smooth_1overP->push_back(m_smoothedBremFit->Get1overP(smoothgraph_fill));
- m_smooth_1overPerr->push_back(m_smoothedBremFit->Get1overPerror(smoothgraph_fill));
- m_smooth_value->push_back(m_smoothedBremFit->GetValue(smoothgraph_fill));
- }
-
- m_forwardparameter_constant = m_forwardparameters.constant;
- m_smoothparameter_constant = m_smoothedparameters.constant;
-
- //Calculate the SeparationQuality
-
- SeparationQuality(m_smoothed_kink,m_forward_kink);
-
- if (m_validationTree && !m_Z_mode)
- m_validationTree->Fill();
- else if (m_validationTree2 && m_Z_mode)
- m_validationTree2->Fill();
-
- delete m_brem_value;
- delete m_brem_phi;
- delete m_brem_theta;
- delete m_brem_energy;
- delete m_brem_UpperBound;
- delete m_brem_LowerBound;
- delete m_forward_kink;
- delete m_smoothed_kink;
- delete m_brem_significance;
- delete m_brem_valueCalibrated;
- delete m_surfaceX;
- delete m_surfaceY;
- delete m_surfaceZ;
-
- delete m_forwardBremFit;
- delete m_smoothedBremFit;
-
-
- delete m_forwardparameter_coefficient;
- delete m_forwardparameter_value;
- delete m_forwardparameter_width;
- delete m_smoothparameter_coefficient;
- delete m_smoothparameter_value;
- delete m_smoothparameter_width;
- delete m_forward_1overP;
- delete m_forward_1overPerr;
- delete m_forward_value;
- delete m_smooth_1overP;
- delete m_smooth_1overPerr;
- delete m_smooth_value;
- delete m_KinkSeparationScores;
- delete m_KinkSeparationScoresErr;
-
-
-}
-
-void Trk::BremFind::CombineParameters()
-{
- //Firstly clear all previous entries
- m_combinedparameters.coefficient.clear();
- m_combinedparameters.width.clear();
- m_combinedparameters.value.clear();
-
- //The combined graph is the backward fit - forward fit
- m_combinedparameters.constant = m_smoothedparameters.constant - m_forwardparameters.constant;
-
- for (int forward_counter(0); forward_counter < (int) m_forwardparameters.value.size(); forward_counter++) {
- m_combinedparameters.coefficient.push_back(-m_forwardparameters.coefficient[forward_counter]);
- m_combinedparameters.value.push_back(m_forwardparameters.value[forward_counter]);
- m_combinedparameters.width.push_back(m_forwardparameters.width[forward_counter]);
- }
-
- for (int smooth_counter(0); smooth_counter < (int) m_smoothedparameters.value.size(); smooth_counter++) {
- m_combinedparameters.coefficient.push_back(m_smoothedparameters.coefficient[smooth_counter]);
- m_combinedparameters.value.push_back(m_smoothedparameters.value[smooth_counter]);
- m_combinedparameters.width.push_back(m_smoothedparameters.width[smooth_counter]);
- }
-
-}
-
-
-//This function returns the value of the tanh series at position x
-double Trk::BremFind::Graph_value(GraphParameter graphparameter, double x)
-{
- double arg = graphparameter.constant;
-
- for (int arg_counter(0); arg_counter < (int) graphparameter.value.size(); arg_counter++) {
- arg+=graphparameter.coefficient[arg_counter]*TMath::TanH(x/graphparameter.width[arg_counter] - graphparameter.value[arg_counter]/graphparameter.width[arg_counter]);
- }
- return arg;
-
-}
-
-//Returns the x value of the local maximum between the min and max terms
-double Trk::BremFind::GetMaximumX(GraphParameter graphparameter, double min, double max)
-{
- double value(0);
- double maximum(0);
- double eval(0);
- for (double value_counter(min); value_counter < max; value_counter++) {
- eval=Graph_value(graphparameter, value_counter);
-
- if (eval > maximum) {
- maximum=eval;
- value=value_counter;
- }
-
- }
- return value;
-}
-
-//Returns the y value of the local maximum between the min and max terms
-double Trk::BremFind::GetMaximumY(GraphParameter graphparameter, double min, double max)
-{
- double maximum(0);
- double eval(0);
- for (double value_counter(min); value_counter < max; value_counter++) {
- eval=Graph_value(graphparameter, value_counter);
- if (eval > maximum) {
- maximum = eval;
- }
- }
- return maximum;
-}
-
-
-void Trk::BremFind::GetBremLocations()
-{
- bool double_back_boolean=false;
- //Normalisation factor introduced because initial parameters were chosen and optimised
- //for 10GeV Pt electrons. So the normalization factor for 10GeV events will be 1.0
-
- double P_T(sin(m_perigee_Theta)*(1/(m_perigee_1overP)));
-
- double normalisation(P_T/10.0);
-
- for (int smooth_counter(0); smooth_counter < (int) m_smoothedparameters.value.size(); smooth_counter++) {
- double_back_boolean=false;
-
- for (int forward_counter(0); forward_counter < (int) m_forwardparameters.value.size(); forward_counter++) {
-
- //backward kink must have a value lower than the forward kink
- if (m_forwardparameters.value[forward_counter] < m_smoothedparameters.value[smooth_counter]) {
- continue;
- }
-
- //If there are multiple back kinks before the forward kink, then the back kinks must
- //iterate forward until it's just before the forward kink
- //Therefore the forward loop is broken to iterate the smooth_counter
- if ((smooth_counter+1 < (int) m_smoothedparameters.value.size()) ) {
- if (m_forwardparameters.value[forward_counter] > m_smoothedparameters.value[smooth_counter+1] && m_smoothedparameters.coefficient[smooth_counter+1]/sin(m_perigee_Theta) > (MIN_KINK_THRESHOLD/normalisation)) {
- double_back_boolean=true;
- }
- }
-
- if (double_back_boolean) {
- break;
- }
-
- //If forward kink is just in front of the backward kink then proceed to find the brem value
- else if (m_forwardparameters.value[forward_counter] > m_smoothedparameters.value[smooth_counter]) {
-
- double brem_separation( GetMaximumY(m_combinedparameters, m_smoothedparameters.value[smooth_counter], m_forwardparameters.value[forward_counter])/sin(m_perigee_Theta) );
-
-
- //Calculates the energy of the brem based on the change in momentum in the
- //1overP parameter
- double coefficient_average( (m_forwardparameters.coefficient[forward_counter]+m_smoothedparameters.coefficient[smooth_counter])/(2*m_perigee_Theta) );
- double p_before_bremF( Graph_value(m_forwardparameters, m_forwardparameters.value[forward_counter] - 2*m_forwardparameters.width[forward_counter])/m_perigee_Theta );
- double p_before_bremS( Graph_value(m_smoothedparameters, m_smoothedparameters.value[smooth_counter] - 2*m_smoothedparameters.width[smooth_counter])/m_perigee_Theta );
- double p_before_bremAverage( (p_before_bremF + p_before_bremS)/2 );
- double p_after_brem( p_before_bremAverage + 2*coefficient_average );
- double brem_amplitude( 1/p_before_bremAverage - 1/p_after_brem );
-
-
- //Check the strength of the brem and its coefficients
- if (brem_amplitude > (BREM_AMPLITUDE*normalisation) && brem_separation > (BREM_SEPARATION/normalisation) && m_forwardparameters.coefficient[forward_counter]/sin(m_perigee_Theta) > (COEFFICIENT_THRESHOLD/normalisation) && m_smoothedparameters.coefficient[smooth_counter]/sin(m_perigee_Theta) > (COEFFICIENT_THRESHOLD/normalisation)) {
-
- m_brem_energy->push_back(brem_amplitude);
- m_brem_value->push_back( GetMaximumX(m_combinedparameters,m_smoothedparameters.value[smooth_counter],m_forwardparameters.value[forward_counter]) );
- m_forward_kink->push_back(m_forwardparameters.value[forward_counter]);
- m_smoothed_kink->push_back(m_smoothedparameters.value[smooth_counter]);
-
- break;
- }
-
- }
-
- }
-
- }
-
-}
-
-//This function obtains the upper and lower bounds of the brem value by scanning the peak
-//till there's a drop in magnitude equal to the desired threshold. If type=1 then, this
-//function will return an upperbound. If type=-1, then this function will return the lower
-//bound. Threshold must be used as a percentage in drop from the maximum Y value.
-double Trk::BremFind::GetBound(double start, int type)
-{
- double maximum(Graph_value(m_combinedparameters,start));
- double eval(maximum);
- double value(start);
-
- while (eval > ((100.0-BOUND_THRESHOLD)/100.0)*maximum) {
- eval = Graph_value(m_combinedparameters,value);
- value+=type;
-
- //if the value is above the end of the inner detector
- if (value > m_smoothedBremFit->LastFitMeasurement()) {
- return m_smoothedBremFit->LastFitMeasurement();
- }
- //if the value is below zero, then return 0 as the lower bound
- if (value < 0.0) {
- return 0.0;
- }
- }
-
- return value;
-}
-
-int Trk::BremFind::ClosestMeasurement(double x, QoverPBremFit *BremFit)
-{
- double min_distance(1000);
- double distance(1000);
- //double value(x);
- int measurement_no(0);
-
- for (int counter(0); counter < BremFit->GetNMeasurements(); counter++) {
- distance = fabs(x - BremFit->GetValue(counter));
- if (distance < min_distance) {
- min_distance = distance;
- //value = BremFit->GetValue(counter);
- measurement_no = counter;
- }
- }
- return measurement_no;
-}
-
-double Trk::BremFind::BremSignificance(double x,QoverPBremFit *forwardBremFit, QoverPBremFit *smoothedBremFit)
-{
- int measurement_noF( ClosestMeasurement(x, forwardBremFit) );
- double PiF( 1/forwardBremFit->Get1overP(measurement_noF) );
-
- //Making sure the next measurement is not on the same surface
- int additional_forward(1);
- double value1F( forwardBremFit->GetValue(measurement_noF+1) );
- while (fabs(x-value1F) < 20.0) {
- additional_forward++;
- value1F = forwardBremFit->GetValue(measurement_noF + additional_forward);
- }
- double Pi1F( 1/forwardBremFit->Get1overP(measurement_noF+1) );
-
-
- int measurement_noS( ClosestMeasurement(x, smoothedBremFit) );
- double PiS( 1/smoothedBremFit->Get1overP(measurement_noS) );
-
- int additional_smoothed(1);
- double value1S( smoothedBremFit->GetValue(measurement_noS-1) );
- while (fabs(x-value1S) < 20.0){
- additional_smoothed++;
- value1S = smoothedBremFit->GetValue(measurement_noS - additional_smoothed);
- }
- double Pi1S( 1/smoothedBremFit->Get1overP(measurement_noS-1) );
-
- double sigf( (1-(Pi1F/PiF))/forwardBremFit->Get1overPerror(measurement_noF) );
- double sigb( (1-(PiS/Pi1S))/smoothedBremFit->Get1overPerror(measurement_noS) );
-
- double sigsum( sigf + sigb);
-
- return std::max(sigsum, std::max(sigf,sigb));
-}
-
-//This function provides an estimate of the phi position using the value
-//The smoothed trajectory is used because it is more accurate in the pixels
-//where the majority of brems occur
-double Trk::BremFind::GetPhiFromValue(double value, QoverPBremFit* smoothedBremFit)
-{
- int closest_measurement_no( ClosestMeasurement(value, smoothedBremFit) );
-
- double value_difference;
- double phi_difference;
-
- if (closest_measurement_no == 0) {
- value_difference = smoothedBremFit->GetValue(closest_measurement_no) - smoothedBremFit->GetValue(closest_measurement_no + 1);
- phi_difference = smoothedBremFit->GetPhi(closest_measurement_no) - smoothedBremFit->GetPhi(closest_measurement_no + 1);
- }
-
- else {
- value_difference = smoothedBremFit->GetValue(closest_measurement_no) - smoothedBremFit->GetValue(closest_measurement_no - 1);
- phi_difference = smoothedBremFit->GetPhi(closest_measurement_no) - smoothedBremFit->GetPhi(closest_measurement_no - 1);
- }
-
- double gradient( phi_difference/value_difference );
- double estimated_phi;
- if (closest_measurement_no == 0)
- estimated_phi = smoothedBremFit->GetPhi(closest_measurement_no+1) + gradient*(value - smoothedBremFit->GetValue(closest_measurement_no + 1));
- else
- estimated_phi = smoothedBremFit->GetPhi(closest_measurement_no-1) + gradient*(value - smoothedBremFit->GetValue(closest_measurement_no - 1));
-
- return estimated_phi;
-
-}
-
-void Trk::BremFind::CalibrateValue()
-{
- double new_value(0);
- double max_value_measurement(m_smoothedBremFit->LastFitMeasurement());
- double min_value_measurement(m_smoothedBremFit->FirstFitMeasurement());
-
-
- for (int counter(0); counter < (int) m_brem_value->size(); counter++) {
-
- double smoothed_kink((*m_smoothed_kink)[counter]);
- double forward_kink((*m_forward_kink)[counter]);
- new_value = ((max_value_measurement - min_value_measurement)-((*m_brem_value)[counter]-min_value_measurement))*smoothed_kink/(max_value_measurement-min_value_measurement);
- new_value += ((*m_brem_value)[counter]-min_value_measurement)*forward_kink/(max_value_measurement-min_value_measurement);
- m_brem_valueCalibrated->push_back( new_value );
- }
-}
-
-
-//This function uses the closest measurement trackparameter and propagates to the brem surface
-//and returns the global position at the intersection. It also fills a vector of trackparameters
-//at the brem surface.
-Amg::Vector3D Trk::BremFind::SurfacePosition(const Trk::TrackParameters& trackparameter, const Trk::Surface& surface, double bremvalue, double bremphi, double bremtheta)
-{
- const Trk::TrackParameters* surfaceParameters1;
- const Trk::TrackParameters* surfaceParameters2;
-
- Amg::Vector3D Pos1 (0, 0, 0);
- Amg::Vector3D Pos2 (0, 0, 0);
- double phidifference1(0.), phidifference2(0.);
-
- //Pushing back just the global position of the brem not the extrapolated
- double x,y,z,r;
- if (m_Z_mode) {
- x = bremvalue * tan(bremtheta) * cos(bremphi);
- y = bremvalue * tan(bremtheta) * sin(bremphi);
- z = bremvalue;
- r = sqrt(x*x+y*y);
- }
- else {
- x = bremvalue * cos(bremphi);
- y = bremvalue * sin(bremphi);
- z = bremvalue * cos(bremtheta);
- r = bremvalue;
- }
-
- const Amg::Vector3D position(x,y,z);
- const Amg::Vector3D momentum((&trackparameter)->momentum());
-
- ATH_MSG_DEBUG("x,y,z,r = " << x << "," << y << "," << z << "," << r);
- ATH_MSG_DEBUG("momentum = " << momentum.mag());
- ATH_MSG_DEBUG("charge = " << (&trackparameter)->charge());
-
- if (!m_usePropagate && momentum.mag() > 0) {
-
- Amg::Transform3D *trans = new Amg::Transform3D(trackparameter.associatedSurface().transform());
- Trk::CylinderSurface cylinderSurface(trans, r, 5000.0);
- const Trk::TrackParameters *bremParameters = new AtaCylinder(position, momentum, (&trackparameter)->charge(), cylinderSurface);
-
- //Broken code - causes FPEs.
- //Trk::CylinderSurface cylinderSurface(new Amg::Transform3D, r, 5000.0);
-
- m_brem_TrackParameters.push_back(bremParameters);
- return position;
- }
-
- //function will propagate both ways and return the trackparameter that has the closest phi
- //to the measured brem phi angle
- surfaceParameters1 = m_propagator->propagate(trackparameter,surface,Trk::oppositeMomentum,true,m_fieldProperties,electron);
- if (surfaceParameters1) {
- Pos1 = surfaceParameters1->position();
- double phi1(atan(Pos1.y()/Pos1.x()));
- phidifference1 = UniqueAngle(phi1 - bremphi);
- }
- surfaceParameters2 = m_propagator->propagate(trackparameter,surface,Trk::alongMomentum,true,m_fieldProperties,electron);
- if (surfaceParameters2) {
- Pos2 = surfaceParameters2->position();
- double phi2((Pos2.y()/Pos2.x()));
- phidifference2 = UniqueAngle(phi2 - bremphi);
- }
- if (surfaceParameters1 && !surfaceParameters2) {
- m_brem_TrackParameters.push_back(surfaceParameters1);
- return Pos1;
- }
- else if (!surfaceParameters1 && surfaceParameters2) {
- m_brem_TrackParameters.push_back(surfaceParameters2);
- return Pos2;
- }
- else if (surfaceParameters1 && surfaceParameters2) {
- if (phidifference1 > phidifference2) {
- m_brem_TrackParameters.push_back(surfaceParameters2);
- return Pos2;
- }
- else {
- m_brem_TrackParameters.push_back(surfaceParameters1);
- return Pos1;
- }
-
- }
-
- //backup method. This really shouldn't be called
- const Trk::TrackParameters* surfaceParameters3 = m_propagator->propagate(trackparameter,surface,Trk::anyDirection,true,m_fieldProperties,electron);
- if (surfaceParameters3) {
- m_brem_TrackParameters.push_back(surfaceParameters3);
- return surfaceParameters3->position();
- }
- m_brem_TrackParameters.push_back(0);
- Amg::Vector3D nullPos(0.0,0.0,0.0);
- return nullPos;
-}
-
-const Trk::Surface* Trk::BremFind::ClosestSurface(double, double, double, const Trk::TrackParameters& trackparameter)
-{
- const Trk::Layer* associatedLayer = m_trackingGeometry->closestMaterialLayer<const Trk::TrackParameters>(trackparameter).object;
- if (associatedLayer) {
- m_brem_layers.push_back(associatedLayer);
- const Trk::Surface* associatedSurface = &(associatedLayer->surfaceRepresentation());
- return associatedSurface;
- }
- else {
- m_brem_layers.push_back(0);
- }
-
- return 0;
-}
-
-
-StatusCode Trk::BremFind::retrieveTrackingGeometry()
-{
- // Retrieve the TrackingGeometry from the DetectorStore
- StatusCode sc = detStore()->retrieve(m_trackingGeometry, m_trackingGeometryName);
- if (sc.isFailure()){
- msg(MSG::FATAL) << "Could not retrieve TrackingGeometry from DetectorStore!" << endmsg;
- }
- return sc;
-
-}
-
-double Trk::BremFind::UniqueAngle(double phi)
-{
- if (phi > TMath::Pi()) {
- return UniqueAngle(phi - TMath::Pi());
- }
- else if (phi < 0.0) {
- return UniqueAngle(phi + TMath::Pi());
- }
- return phi;
-}
-
-int Trk::BremFind::GetNBrems()
-{
- return m_nBrems;
-}
-
-void Trk::BremFind::MakeTrackStateOnSurfaces ()
-{
-
- for (int make_counter(0); make_counter < m_nBrems; make_counter++) {
-
- std::bitset<TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
- typePattern.set(TrackStateOnSurface::BremPoint);
- const Trk::TrackStateOnSurface* BremTrackStateOnSurface = new Trk::TrackStateOnSurface(0,m_brem_TrackParameters[make_counter],0,m_EstimatedBremOnTrack[make_counter],typePattern);
-
- if (!m_Z_mode)
- m_brem_trackStateOnSurface.push_back(BremTrackStateOnSurface);
- else
- delete BremTrackStateOnSurface;
-
- }
-
-}
-
-const Trk::TrackStateOnSurface* Trk::BremFind::GetEstimatedBremOnTrack(int n)
-{
- return m_brem_trackStateOnSurface[n];
-}
-
-
-//A function that rates the quality of a bremsstrahlung signal based on the separation quality
-//of the forward and smoothed fit. A lower score indicates that the measurements are unlikely to
-//be the same value. Therefore a lower score means a better fit.
-
-void Trk::BremFind::SeparationQuality(std::vector<double> *backwardkink, std::vector<double> *forwardkink)
-{
- int ForwardHit_Counter(0);
- int SmoothHit_Counter(1);
- std::vector<double> OverlapScore(0);
- std::vector<double> OverlapErr(0);
-
-
- for (int kink_counter(0); kink_counter < (int) backwardkink->size(); kink_counter++) {
-
- OverlapScore.clear();
- OverlapErr.clear();
-
- ForwardHit_Counter = 0;
- SmoothHit_Counter = 1;
-
- for (; ForwardHit_Counter < (int) m_forward_1overP->size(); ForwardHit_Counter++,SmoothHit_Counter++) {
-
-
- if ( (*m_forward_value)[ForwardHit_Counter] > (*backwardkink)[kink_counter] && (*m_forward_value)[ForwardHit_Counter] < (*forwardkink)[kink_counter] && (*m_smooth_value)[SmoothHit_Counter] > (*backwardkink)[kink_counter] && (*m_smooth_value)[SmoothHit_Counter] < (*forwardkink)[kink_counter]) {
-
- std::pair<double,double> probabilitypair(ProbabilityScore((*m_forward_1overP)[ForwardHit_Counter],3*(*m_forward_1overPerr)[ForwardHit_Counter],(*m_smooth_1overP)[SmoothHit_Counter],3*(*m_smooth_1overPerr)[SmoothHit_Counter]));
- OverlapScore.push_back(probabilitypair.first);
- OverlapErr.push_back(probabilitypair.second);
-
- }
-
- }
-
- //Now combing the multiple measurements into one score
-
- double sum1(0);
- double sum2(0);
- double LowestScore(1.0);
- double LowestScoreErr(0.0);
-
- for(int score_counter(0); score_counter < (int) OverlapScore.size(); score_counter++) {
-
- sum1 += OverlapScore[score_counter]/(pow(OverlapErr[score_counter],2));
- sum2 += 1/(pow(OverlapErr[score_counter],2));
-
- if (OverlapScore[score_counter] < LowestScore || score_counter==0) {
- LowestScore = OverlapScore[score_counter];
- LowestScoreErr = OverlapErr[score_counter];
- }
- }
-
- if (sum2==0) {
- m_KinkSeparationScores->push_back(1.0);
- m_KinkSeparationScoresErr->push_back(0.0);
- }
- else {
- m_KinkSeparationScores->push_back(sum1/sum2);
- m_KinkSeparationScoresErr->push_back(LowestScoreErr);
- //m_KinkSeparationScoresErr->push_back(1/sqrt(sum2));
- }
-
- }
-
-}
-
-
-std::pair<double,double> Trk::BremFind::ProbabilityScore(double Gauss1_mean, double Gauss1_var, double Gauss2_mean, double Gauss2_var)
-{
- double value(exp(-(Gauss1_mean - Gauss2_mean)*(Gauss1_mean - Gauss2_mean)/(2*(Gauss1_var*Gauss1_var + Gauss2_var*Gauss2_var)))*(1/Gauss1_var*Gauss2_var*2*sqrt(TMath::Pi()))*sqrt(2*Gauss1_var*Gauss1_var*Gauss2_var*Gauss2_var/(Gauss1_var*Gauss1_var + Gauss2_var*Gauss2_var)));
- double error(value*(Gauss1_mean - Gauss2_mean)*(Gauss1_var + Gauss2_var)/(Gauss1_var*Gauss1_var + Gauss2_var*Gauss2_var));
-
- std::pair<double,double> probabilityscore;
- probabilityscore.first = value;
- probabilityscore.second = fabs(error);
-
- return probabilityscore;
-
-}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
index 0bebf8563fe0a42c2410eb4218b629239a54be40..b7c020558af810a90fa10d0a96760e36e67955d5 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
@@ -18,7 +18,6 @@ decription : Implementation code for Gaussian Sum Fitter class
#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "TrkGaussianSumFilter/IForwardGsfFitter.h"
#include "TrkGaussianSumFilter/IGsfSmoother.h"
-#include "TrkGaussianSumFilter/IBremsstrahlungFinder.h"
#include "TrkEventUtils/MeasurementBaseComparisonFunction.h"
#include "TrkEventUtils/PrepRawDataComparisonFunction.h"
@@ -61,13 +60,10 @@ Trk::GaussianSumFitter::GaussianSumFitter(const std::string& type, const std::st
m_reintegrateOutliers(false),
m_makePerigee(true),
m_directionToPerigee(Trk::oppositeMomentum),
- m_runBremFinder(false),
m_refitOnMeasurementBase(true),
m_doHitSorting(true),
m_trkParametersComparisonFunction(0),
m_stateCombiner("Trk::MultiComponentStateCombiner"),
- m_BremFind("Trk::BremFind"),
- m_BremFind2("Trk::BremFind"),
m_chronoSvc("ChronoStatSvc", name),
m_inputPreparator(0),
m_FitPRD(0),
@@ -130,9 +126,6 @@ Trk::GaussianSumFitter::GaussianSumFitter(const std::string& type, const std::st
declareProperty("SortingReferencePoint", m_sortingReferencePoint );
declareProperty("StateCombiner", m_stateCombiner );
declareProperty("ValidationMode", m_validationMode );
- declareProperty("BremFind", m_BremFind );
- declareProperty("BremFind2", m_BremFind2 );
- declareProperty("runBremFinder", m_runBremFinder );
declareProperty("GsfSmoother",m_gsfSmoother);
declareProperty("ForwardGsfFitter",m_forwardGsfFitter);
declareProperty("EventInfoKey", m_readKey="EventInfo");
@@ -169,17 +162,6 @@ StatusCode Trk::GaussianSumFitter::initialize()
// Request the state combiner
ATH_CHECK ( m_stateCombiner.retrieve() );
- //Request the brem finder
- if (m_runBremFinder){
- if ( m_BremFind.retrieve().isFailure() || m_BremFind2.retrieve().isFailure() ) {
- msg(MSG::WARNING) << "Request is to retrieve the bremsstrahlung finder failed... turning off brem finding!" << endmsg;
- m_runBremFinder = false;
- }
- }
- else {
- m_BremFind.disable();
- m_BremFind2.disable();
- }
// Request the RIO_OnTrack creator
// No need to return if RioOnTrack creator tool, only if PrepRawData is used in fit
@@ -573,16 +555,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
}
}
- //Find bremsstrahlung points and save them
- if (m_runBremFinder) {
- if (msgLvl(MSG::DEBUG)) msg() << "Entering the BremFind tool" << endmsg;
-
- //The z mode must be ahead of the r mode to make sure the right TSOS is being pushed into the smoothedTrajectory
- m_BremFind2->BremFinder(*forwardTrajectory, *smoothedTrajectory,true);
- m_BremFind->BremFinder(*forwardTrajectory, *smoothedTrajectory,false);
-
- }
-
+
// Delete forward trajectory. New memory was assigned in ForwardGsfFitter.
delete forwardTrajectory;
@@ -594,15 +567,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::PrepRawDataSet& prepRawD
//Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
- //Add brem information to Track
- if (m_runBremFinder) {
- //Filling in the brem information in the smoothedTrajectory
- for (int brem_counter(0); brem_counter < m_BremFind->GetNBrems(); brem_counter++) {
- smoothedTrajectory->push_back(m_BremFind->GetEstimatedBremOnTrack(brem_counter));
- }
- }
-
-
+
// Create new track
Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
info.setTrackProperties(TrackInfo::BremFit);
@@ -753,15 +718,6 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
}
}
- //Find bremsstrahlung points and save them
- if (m_runBremFinder) {
- if (msgLvl(MSG::DEBUG)) msg() << "Entering the BremFind tool" << endmsg;
-
- //The z mode must be ahead of the r mode to make sure the right TSOS is being pushed into the smoothedTrajectory
- m_BremFind2->BremFinder(*forwardTrajectory, *smoothedTrajectory,true);
- m_BremFind->BremFinder(*forwardTrajectory, *smoothedTrajectory,false);
- }
-
//Delete forward trajectory. New memory was assigned in ForwardGsfFitter.
delete forwardTrajectory;
@@ -774,14 +730,7 @@ Trk::Track* Trk::GaussianSumFitter::fit ( const Trk::MeasurementSet& measurem
//Reverse the order of the TSOS's to make be order flow from inside to out
std::reverse(smoothedTrajectory->begin(), smoothedTrajectory->end());
- //Add brem information to the track
- if (m_runBremFinder) {
- //Filling in the brem information in the smoothedTrajectory
- for (int brem_counter(0); brem_counter < m_BremFind->GetNBrems(); brem_counter++) {
- smoothedTrajectory->push_back(m_BremFind->GetEstimatedBremOnTrack(brem_counter));
- }
- }
-
+
// Create new track
Trk::TrackInfo info(Trk::TrackInfo::GaussianSumFilter, particleHypothesis);
info.setTrackProperties(TrackInfo::BremFit);
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.cxx
deleted file mode 100644
index af6df4cb5a43b0db0cde733317a51dd37f47b7dd..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.cxx
+++ /dev/null
@@ -1,711 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/********************************************************************************
- QoverPBremFit.cxx * description
-
-completed : Tuesday 18th November 2008
-author : Tony Shao
-email : Qi.Tao.Shao@cern.ch
-description : Class called by BremFind to fit 1overp vs R plots using a sum of Tanh graphs
-**********************************************************************************/
-#include "QoverPBremFit.h"
-
-#include "TrkGaussianSumFilter/IMultiComponentStateCombiner.h"
-#include "TrkMultiComponentStateOnSurface/MultiComponentStateOnSurface.h"
-#include "TrkParameters/TrackParameters.h"
-#include "TrkTrack/Track.h"
-#include "TrkFitterUtils/FitterTypes.h"
-
-#include <math.h>
-#include <algorithm>
-#include "TMath.h"
-
-
-
-Trk::QoverPBremFit::QoverPBremFit(int type, bool Z_mode):
- m_stateCombiner("Trk::MultiComponentStateCombiner")
-{
- m_potentialvalue = new std::vector<double>;
- m_potentialwidth = new std::vector<double>;
- m_1overPvalue = new std::vector<double>;
- m_1overPvalueerror = new std::vector<double>;
- m_PhiValue = new std::vector<double>;
- m_PhiValueerror = new std::vector<double>;
- m_ThetaValue = new std::vector<double>;
- m_ThetaValueerror = new std::vector<double>;
- m_RadiusValue = new std::vector<double>;
- m_Zvalue = new std::vector<double>;
- m_FitValue = new std::vector<double>;
- m_trackParameters = new std::vector<const Trk::TrackParameters*>;
- m_type = type;
- m_Z_mode = Z_mode;
- m_surfaceCounter = 0;
- m_charge =0;
-}
-
-StatusCode Trk::QoverPBremFit::initialize()
-{
- // Request the state combiner
- if ( m_stateCombiner.retrieve().isFailure() ){
- printf("Failed to retrieve m_stateCombiner");
- return StatusCode::FAILURE;
- }
- return StatusCode::SUCCESS;
-}
-
-
-Trk::QoverPBremFit::~QoverPBremFit()
-{
- for (unsigned int d(0);d < m_trackParameters->size(); d++) {
- delete (*m_trackParameters)[d];
- }
-
- delete m_potentialvalue;
- delete m_potentialwidth;
- delete m_1overPvalue;
- delete m_1overPvalueerror;
- delete m_PhiValue;
- delete m_PhiValueerror;
- delete m_ThetaValue;
- delete m_ThetaValueerror;
- delete m_RadiusValue;
- delete m_Zvalue;
- delete m_FitValue;
- delete m_trackParameters;
-}
-
-
-
-Trk::GraphParameter Trk::QoverPBremFit::GetParameters(const DataVector <const Trk::TrackStateOnSurface>& Trajectory)
-{
- FillVariables(Trajectory);
-
- for (int copy_counter(0); (int) copy_counter < GetNMeasurements(); copy_counter++) {
- if (m_Z_mode)
- m_FitValue->push_back(fabs((*m_Zvalue)[copy_counter]));
- else
- m_FitValue->push_back((*m_RadiusValue)[copy_counter]);
- }
-
-
- FindPotentials();
-
- m_graphparameter=FindParameters();
-
- m_mergedparameter = MergeParameters(m_graphparameter);
-
-
- return m_mergedparameter;
-}
-
-
-void
-Trk::QoverPBremFit::FillVariables(const DataVector <const Trk::TrackStateOnSurface>& Trajectory) {
- //iterators will be constructed as NULL
- DataVector<const Trk::TrackStateOnSurface>::const_iterator trackStateOnSurface;
- DataVector<const Trk::TrackStateOnSurface>::const_iterator Trajectory_end;
-
-
- if (m_type==1) {
- trackStateOnSurface = Trajectory.begin();
- Trajectory_end = Trajectory.end();
- } else if (m_type==-1) {
- trackStateOnSurface = Trajectory.end();
- trackStateOnSurface--;
- Trajectory_end = Trajectory.begin();
- Trajectory_end--;
- } else {
- //trackStateOnSurface and Trajectory_end will be NULL in this case, as constructed
- return;
- }
-
-
- m_surfaceCounter=0;
-
-
- for ( ; trackStateOnSurface != Trajectory_end; ) {
- const Trk::MultiComponentState *TrajectoryMultiState=0;
- const Trk::MultiComponentStateOnSurface *TrajectoryMultiStateOnSurface = dynamic_cast<const Trk::MultiComponentStateOnSurface*>(*trackStateOnSurface);
- if (m_surfaceCounter < TRKFGSF_VALSURFACES) {
- if (!TrajectoryMultiStateOnSurface) {
- // Create new multiComponentState from a single state
- Trk::ComponentParameters componentParameters( (*trackStateOnSurface)->trackParameters(), 1.0);
- TrajectoryMultiState = new Trk::MultiComponentState(componentParameters);
- } else {
- TrajectoryMultiState = TrajectoryMultiStateOnSurface->components();
- }
-
- const Amg::Vector3D posOnSurf = TrajectoryMultiStateOnSurface->trackParameters()->position();
-
-
- //Fill the radius values
- m_RadiusValue->push_back( posOnSurf.perp() );
- m_Zvalue->push_back( posOnSurf.z() );
-
- //Specify in joboptions whether to use mode or mean
- const Trk::TrackParameters* CombinedMultiState = m_stateCombiner->combine(*TrajectoryMultiState,true);
-
- if (CombinedMultiState->parameters()[Trk::qOverP] > 0.0) m_charge = 1.0;
- else m_charge = -1.0;
-
- //Fill the 1overPvalues in GeV
- m_trackParameters->push_back( CombinedMultiState );
- m_1overPvalue->push_back( 1000.0*fabs(CombinedMultiState->parameters()[Trk::qOverP]) );
- m_PhiValue->push_back( CombinedMultiState->parameters()[Trk::phi] );
- m_ThetaValue->push_back( CombinedMultiState->parameters()[Trk::theta] );
-
- const AmgSymMatrix(5)* measuredCov = CombinedMultiState->covariance();
-
- if (measuredCov) {
- m_1overPvalueerror->push_back( 1000.0*sqrt((*measuredCov)(Trk::qOverP,Trk::qOverP)));
- m_PhiValueerror->push_back( sqrt((*measuredCov)(Trk::phi,Trk::phi)));
- m_ThetaValueerror->push_back( sqrt((*measuredCov)(Trk::theta,Trk::theta)));
- }
- else {
- m_1overPvalueerror->push_back( 1000.0*fabs(CombinedMultiState->parameters()[Trk::qOverP]) );
- m_PhiValueerror->push_back( fabs(CombinedMultiState->parameters()[Trk::phi]) );
- m_ThetaValueerror->push_back( fabs(CombinedMultiState->parameters()[Trk::theta]) );
- }
- }
-
- if (!TrajectoryMultiStateOnSurface) {
- delete TrajectoryMultiState;
- }
- m_surfaceCounter++;
- if (m_type==1) {
- trackStateOnSurface++;
- }
- else if (m_type==-1) {
- trackStateOnSurface--;
- }
- }//end of for loop
-}
-
-
-void Trk::QoverPBremFit::FindPotentials()
-{
- double width;
- double oneoverPgradient;
-
- //Search for potential brem locations
- for (int datapoint=0; datapoint < m_surfaceCounter-1; datapoint++) {
- //Calculates gradient between adajacent points
- oneoverPgradient=(1 - (*m_1overPvalue)[datapoint]/(*m_1overPvalue)[datapoint+1])/((*m_FitValue)[datapoint+1]-(*m_FitValue)[datapoint]);
-
- if (oneoverPgradient > POTENTIAL_BREM_THRESHOLD) {
-
- m_potentialvalue->push_back(((*m_FitValue)[datapoint+1]+(*m_FitValue)[datapoint])/2);
- width=((*m_FitValue)[datapoint+1]-(*m_FitValue)[datapoint])/2;
-
- if (width < 20) {
- width=20;
- }
-
- m_potentialwidth->push_back(width);
-
- }
- }
-}
-
-
-Trk::GraphParameter Trk::QoverPBremFit::FindParameters()
-{
-
- GraphParameter graphparameter;
- int nparameters = (int) m_potentialvalue->size() + 1;
- int totalparameters = nparameters;
- int nsurfaces = m_1overPvalue->size();
-
-
- std::vector<Element> elements;
- elements.clear();
- Element thePoint;
-
- //Fill in the elements
- for (int elementfill=0; elementfill < totalparameters-1; elementfill++) {
-
- thePoint.value=(*m_potentialvalue)[elementfill];
- thePoint.width=(*m_potentialwidth)[elementfill];
- thePoint.sign=true;
- elements.push_back(thePoint);
- }
-
- while (true) {
- //printf("Starting the matrix iterations\n");
- Amg::VectorX theta(nparameters); theta.setZero();
- Amg::VectorX ymatrix(nsurfaces); ymatrix.setZero();
- Amg::MatrixX Amatrix (nsurfaces,nparameters); Amatrix.setZero();
- Amg::MatrixX Gy (nsurfaces,nsurfaces); Gy.setZero();
-
- //temporary matrices to store the operations
- Amg::MatrixX ATranspose(nparameters,nsurfaces); ATranspose.setZero();
- Amg::MatrixX AGA(nparameters,nparameters); AGA.setZero();
- Amg::MatrixX AGAinverse(nparameters,nparameters); AGAinverse.setZero();
-
- //printf("nparameters=%i\n",nparameters);
-
-
- //Filling each of the matrices --- y matrix first
- for (int yfill=0; yfill < nsurfaces; yfill++) {
- ymatrix[yfill]=(*m_1overPvalue)[yfill];
- }
-
-
- //Filling the A matrix.
- for (int afill=0; afill < nsurfaces; afill++) {
-
- int afill3(0);
-
- //afill2 starts at -1 because we want the elements vector to start at zero.
- //The first iteration is taken away by the constant parameters.
- for (int afill2=-1; afill2 < totalparameters-1; afill2++) {
-
- //First column of the A matrix is a constant
- if (afill3==0){
- Amatrix(afill,afill3)=1;
- afill3++;
- }
-
- //Only the terms with positive coefficients are kept
- else if (elements[afill2].sign) {
-
- Amatrix(afill,afill3)=TMath::TanH((*m_FitValue)[afill]/elements[afill2].width-elements[afill2].value/elements[afill2].width);
- afill3++;
- }
-
- }
- }
-
-
- //Fill the Gy matrix (the error matrix)
- for (int gfill(0); gfill < nsurfaces; gfill++) {
-
- for (int gfill2(0); gfill2 < nsurfaces; gfill2++) {
-
- if (gfill==gfill2) {
- Gy(gfill,gfill2)=1/(*m_1overPvalueerror)[gfill];
- }
- else {
- Gy(gfill,gfill2)=0;
- }
- }
- }
-
-
-
- //int inv;
- ATranspose=Amatrix.transpose();
- AGA=ATranspose*(Gy*(Amatrix));
- AGAinverse=AGA.inverse();
- //printf("Inv=%i\n",inv);
-
- theta=AGAinverse*(ATranspose*(Gy*ymatrix));
-
-
- //Check the signs of the calculated coefficients. Negative coefficients will be taken away
- //on the next iteration until all coefficients are >= 0
-
- int nNegative(0);
-
- int signcheck2(1); //Iterations start at 1 because the constant parameter doesn't matter
- for (int signcheck(0); signcheck < totalparameters-1; signcheck++) {
-
- if (!elements[signcheck].sign) { //if the term is already taken away, then ignore it
- continue;
- }
- else if (theta[signcheck2] < 0) { //checks to the sign of the calculated coefficient
- elements[signcheck].sign=false;
- nNegative++; //counts the number of negative coefficients
- signcheck2++;
- }
- else {
- signcheck2++;
- }
- }
-
- if (nNegative==0) {//Ends the iterations because all terms are positive
-
- int graphfill2(0);
- for(int graphfill(-1); graphfill < totalparameters-1; graphfill++) {
- if (graphfill2==0) { //Fills the constant with the corresponding theta element
- graphparameter.constant=theta[graphfill2];
- if (isnan(graphparameter.constant))
- graphparameter.constant=0.0;
- graphfill2++;
- }
-
- else if (!elements[graphfill].sign) {
- //if the parameter is negative then its coefficient is assigned zero
- graphparameter.value.push_back(elements[graphfill].value);
- graphparameter.coefficient.push_back(0.0);
- graphparameter.width.push_back(20.0);
- }
-
- else {
- //fills the graphparameter with the coefficients, radii and width
- if (isnan(theta[graphfill2])) {
- graphparameter.coefficient.push_back(0.0);
- graphparameter.value.push_back(elements[graphfill].value);
- graphparameter.width.push_back(elements[graphfill].width);
- }
- else {
- graphparameter.value.push_back(elements[graphfill].value);
- graphparameter.coefficient.push_back(theta[graphfill2]);
- graphparameter.width.push_back(elements[graphfill].width);
- }
- graphfill2++;
- }
- }
-
- break;
-
- }
- else {
- nparameters-=nNegative; //resets the matrices with less degrees of freedom
- }
- }
-
- return graphparameter;
-}
-
-
-Trk::GraphParameter Trk::QoverPBremFit::MergeParameters(GraphParameter graphparameter)
-{
- GraphParameter MergedParameters;
- MergedParameters.constant=graphparameter.constant;
-
- double joinscore(0);
- std::vector<int> nonzero_coefficients; //contains the coefficient numbers that are non-zero
-
-
- for (int coeff_test(0); coeff_test < (int) graphparameter.coefficient.size(); coeff_test++) {
- if ( graphparameter.coefficient[coeff_test] != 0) {
- nonzero_coefficients.push_back(coeff_test);
- }
- }
-
- //Returns the singular case
- if ( (int) nonzero_coefficients.size() == 1 ) {
- MergedParameters.coefficient.push_back(graphparameter.coefficient[nonzero_coefficients[0]]);
- MergedParameters.value.push_back(graphparameter.value[nonzero_coefficients[0]]);
- MergedParameters.width.push_back(graphparameter.width[nonzero_coefficients[0]]);
- return MergedParameters;
- }
-
- std::vector<bool> joinboolean ((int) nonzero_coefficients.size()-1,false);
- double kinkdistance; //distance between adjacent kinks
- double maxmeasure; //the maximum distance between measurements between the adjacent kinks
- //int npotential; //number of potential brem points between kinks
- int nmeasurement; //number of measurements between kinks
- double regionbias; //takes into account the joining bias of the region based on value
-
- //The first element of the joinboolean array will represent whether
- //kink #1 and kink#2 should be joined
-
- for (int join_counter(0); join_counter < (int) nonzero_coefficients.size()-1; join_counter++) {
-
- joinscore=0;
-
- kinkdistance = graphparameter.value[nonzero_coefficients[join_counter+1]] - graphparameter.value[nonzero_coefficients[join_counter]];
- maxmeasure = maxvalue(graphparameter.value[nonzero_coefficients[join_counter]],graphparameter.value[nonzero_coefficients[join_counter+1]]);
- //npotential = npotentials(graphparameter.value[nonzero_coefficients[join_counter]],graphparameter.value[nonzero_coefficients[join_counter+1]]);
- nmeasurement = nmeasurements(graphparameter.value[nonzero_coefficients[join_counter]],graphparameter.value[nonzero_coefficients[join_counter+1]]);
- regionbias = (LastFitMeasurement()/2-graphparameter.value[nonzero_coefficients[join_counter+1]]+LastFitMeasurement()/2-graphparameter.value[nonzero_coefficients[join_counter]])/5;
-
- joinscore+= -10*nmeasurement;
- joinscore+= 2*maxmeasure;
- joinscore+= 50-kinkdistance;
- joinscore+= m_type*regionbias;
- //joinscore+= 20*npotential;
-
- //printf("Join Score%i=%f\n",join_counter,joinscore);
-
- //Scoring system thresholds
- if (kinkdistance > 350) {
- joinboolean[join_counter] = false;
- continue;
- }
-
- if (joinscore > 0) {
- joinboolean[join_counter] = true;
- }
-
- else {
- joinboolean[join_counter] = false;
- }
-
- }
-
- //Use the joinboolean array to join multiple smaller kinks into one larger one
-
- int nkinks_merged(0);
-
- std::vector<double> coefficient_vector;
- std::vector<double> width_vector;
- std::vector<double> value_vector;
-
- for (int merge_counter(0); merge_counter < (int) nonzero_coefficients.size(); merge_counter++) {
-
- if (merge_counter < (int) nonzero_coefficients.size()-1) {
-
- if (joinboolean[merge_counter]) {
- coefficient_vector.push_back(graphparameter.coefficient[nonzero_coefficients[merge_counter]]);
- width_vector.push_back(graphparameter.width[nonzero_coefficients[merge_counter]]);
- value_vector.push_back(graphparameter.value[nonzero_coefficients[merge_counter]]);
- nkinks_merged++;
- }
-
- else if (!joinboolean[merge_counter]) {
- coefficient_vector.push_back(graphparameter.coefficient[nonzero_coefficients[merge_counter]]);
- value_vector.push_back(graphparameter.value[nonzero_coefficients[merge_counter]]);
- width_vector.push_back(graphparameter.width[nonzero_coefficients[merge_counter]]);
- MergedParameters.coefficient.push_back(coefficient_sum(coefficient_vector));
- MergedParameters.value.push_back(value_sum(value_vector,coefficient_vector));
- MergedParameters.width.push_back(width_sum(width_vector,coefficient_vector,value_vector));
- coefficient_vector.clear();
- value_vector.clear();
- width_vector.clear();
- nkinks_merged=0;
- }
- }
-
- else {
- coefficient_vector.push_back(graphparameter.coefficient[nonzero_coefficients[merge_counter]]);
- value_vector.push_back(graphparameter.value[nonzero_coefficients[merge_counter]]);
- width_vector.push_back(graphparameter.width[nonzero_coefficients[merge_counter]]);
- MergedParameters.coefficient.push_back(coefficient_sum(coefficient_vector));
- MergedParameters.value.push_back(value_sum(value_vector,coefficient_vector));
- MergedParameters.width.push_back(width_sum(width_vector,coefficient_vector,value_vector));
- }
-
- }
-
- return MergedParameters;
-}
-
-
-//This function calculates the number of measurements between kinks
-int Trk::QoverPBremFit::nmeasurements (double start, double end)
-{
- int npoints(0);
- for (int i=0; i < (int) m_FitValue->size(); i++) {
- if ((*m_FitValue)[i] > start && (*m_FitValue)[i] < end) {
- npoints++;
- }
- }
-
- return npoints;
-}
-
-//This function calculates the number of potential brem points between kinks
-int Trk::QoverPBremFit::npotentials(double start, double end)
-{
- int potentialcount(0);
- for (int i=0; i < (int) m_potentialvalue->size(); i++) {
- if ( (*m_potentialvalue)[i] > start && (*m_potentialvalue)[i] < end) {
- potentialcount++;
- }
-
- }
- return potentialcount;
-}
-
-
-//This function calculates the maximum distance between measurements within adjacent kinks
-//Note: the measured points considered can be outside the radii of the kinks, hence the use
-//of "if" conditions #2 and #4 below.
-double Trk::QoverPBremFit::maxvalue(double start, double end)
-{
- double max(-1);
- double max2(-1);
- int endcounter(0);
- for (int value_counter(0); value_counter < (int) m_FitValue->size(); value_counter++) {
-
- if ((*m_FitValue)[value_counter] < start) {
- continue;
- }
-
- if ((*m_FitValue)[value_counter] > end && endcounter > 0) {
- continue;
- }
-
- //A measurement outside the value of the kink will counter if it is the first and only one
- else if ((*m_FitValue)[value_counter] > end && endcounter==0) {
- max2 = (*m_FitValue)[value_counter] - (*m_FitValue)[value_counter - 1];
- endcounter++;
- }
-
- else {
- max2 = (*m_FitValue)[value_counter] - (*m_FitValue)[value_counter - 1];
- }
-
- if (max2 > max) {
- max = max2;
- }
-
- }
- return max;
-
-}
-
-
-double Trk::QoverPBremFit::coefficient_sum (std::vector<double> coefficients)
-{
- double sum(0);
- for (int sum_counter(0); sum_counter < (int) coefficients.size(); sum_counter++) {
- sum+=coefficients[sum_counter];
- }
-
- return sum;
-
-}
-
-double Trk::QoverPBremFit::value_sum (std::vector<double> radii, std::vector<double> coefficients)
-{
- double value(0);
- double mod_squared(0);
- for (int value_counter(0); value_counter < (int) radii.size(); value_counter++) {
- value+=coefficients[value_counter]*coefficients[value_counter]*radii[value_counter];
- mod_squared+=coefficients[value_counter]*coefficients[value_counter];
- }
- value/=mod_squared;
- return value;
-}
-
-double Trk::QoverPBremFit::width_sum (std::vector<double> widths, std::vector<double> coefficients, std::vector<double> radii)
-{
- double width(0);
- double coefficient_average(0);
- double mod_squared(0);
-
- //singular case
- if ( (int)coefficients.size() == 1)
- {
- return widths[0];
- }
-
-
- for (int coefficient_counter(0); coefficient_counter < (int)coefficients.size(); coefficient_counter++) {
- coefficient_average+=coefficients[coefficient_counter]/coefficients.size();
- mod_squared+=coefficients[coefficient_counter]*coefficients[coefficient_counter];
- }
-
- for (int width_counter(0); width_counter < (int) widths.size()-1 ; width_counter++) {
- width+=std::min(coefficients[width_counter]*coefficients[width_counter]*(radii[width_counter+1]-radii[width_counter]),coefficients[width_counter+1]*coefficients[width_counter+1]*(radii[width_counter+1]-radii[width_counter]));
- }
- width/=mod_squared;
-
- for (int width_counter2(0); width_counter2 < (int) widths.size(); width_counter2++) {
- width+=coefficients[width_counter2]*coefficients[width_counter2]*widths[width_counter2]/(coefficients[width_counter2]*coefficients[width_counter2]+coefficient_average*coefficient_average);
- }
-
- return width;
-}
-
-double Trk::QoverPBremFit::GetPerigee_1overP()
-{
- return (*m_1overPvalue)[0];
-}
-
-double Trk::QoverPBremFit::GetPerigee_Phi()
-{
- return (*m_PhiValue)[0];
-}
-
-double Trk::QoverPBremFit::GetPerigee_Theta()
-{
- return (*m_ThetaValue)[0];
-}
-
-double Trk::QoverPBremFit::GetPerigee_d0()
-{
- return (*m_RadiusValue)[0];
-}
-
-double Trk::QoverPBremFit::GetPerigee_z0()
-{
- return (*m_Zvalue)[0];
-}
-
-double Trk::QoverPBremFit::Get1overP(int n)
-{
- return (*m_1overPvalue)[n];
-}
-
-double Trk::QoverPBremFit::GetPhi(int n)
-{
- return (*m_PhiValue)[n];
-}
-
-double Trk::QoverPBremFit::GetTheta(int n)
-{
- return (*m_ThetaValue)[n];
-}
-
-double Trk::QoverPBremFit::GetRadius(int n)
-{
- return (*m_RadiusValue)[n];
-}
-
-double Trk::QoverPBremFit::GetZ(int n)
-{
- return (*m_Zvalue)[n];
-}
-
-double Trk::QoverPBremFit::Get1overPerror(int n)
-{
- return (*m_1overPvalueerror)[n];
-}
-
-double Trk::QoverPBremFit::GetPhierror(int n)
-{
- return (*m_PhiValueerror)[n];
-}
-
-double Trk::QoverPBremFit::GetThetaerror(int n)
-{
- return (*m_ThetaValueerror)[n];
-}
-
-int Trk::QoverPBremFit::GetNMeasurements()
-{
- return (int) (*m_1overPvalue).size();
-}
-
-double Trk::QoverPBremFit::GetCharge()
-{
- return m_charge;
-}
-
-const Trk::TrackParameters* Trk::QoverPBremFit::GetTrackParameters(int n)
-{
- return (*m_trackParameters)[n];
-}
-
-int Trk::QoverPBremFit::GetnSurfaces()
-{
- return m_surfaceCounter;
-}
-
-double Trk::QoverPBremFit::GetValue(int n)
-{
- return (*m_FitValue)[n];
-}
-
-
-double Trk::QoverPBremFit::LastFitMeasurement()
-{
- int lastmeasurement(GetNMeasurements() - 1);
- return (*m_FitValue)[lastmeasurement];
-}
-
-double Trk::QoverPBremFit::FirstFitMeasurement()
-{
- if (m_type == 1)
- return (*m_FitValue)[0];
- else if (m_type == -1)
- return (*m_FitValue)[1];
- return (*m_FitValue)[0];
-}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.h b/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.h
deleted file mode 100644
index 91d187aaf16b3a86b4d0cc4b53213f05123063a9..0000000000000000000000000000000000000000
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/QoverPBremFit.h
+++ /dev/null
@@ -1,129 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/********************************************************************************
- QoverPBremFit.h * description
-
- completed : Tuesday 18th November 2008
- author : Tony Shao
- email : Qi.Tao.Shao@cern.ch
- description : Class called by BremFind to fit QoverP vs R plots using a sum of Tanh graphs
- **********************************************************************************/
-
-#ifndef Trk_QoverPBremFit_H
-#define Trk_QoverPBremFit_H
-
-#include <vector>
-
-#include "GaudiKernel/ToolHandle.h"
-#include "AthContainers/DataVector.h"
-#include "TrkParameters/TrackParameters.h"
-#include "TrkGaussianSumFilter/BremFind.h"
-
-namespace Trk {
-
- class IMultiComponentStateCombiner;
- class MultiComponentStateOnSurface;
- class TrackStateOnSurface;
-
-
-#define TRKFGSF_VALSURFACES 100
-#define TRKGSF_VALSTATES 24
-#define POTENTIAL_BREM_THRESHOLD 0.0001
-
-
- class QoverPBremFit {
-
- // forbid copying
- QoverPBremFit(const QoverPBremFit &) = delete ;
- QoverPBremFit &operator=(const QoverPBremFit &) = delete;
-
- public:
-
- QoverPBremFit(int, bool); //type=1 for forward fit || type=-1 for backward fit
- //bool for Z mode or not. True for Z mode.
-
- StatusCode initialize();
-
- ~QoverPBremFit();
-
- GraphParameter GetParameters(const DataVector <const Trk::TrackStateOnSurface>&);
-
- double GetPerigee_1overP();
- double GetPerigee_Phi();
- double GetPerigee_Theta();
- double GetPerigee_d0();
- double GetPerigee_z0();
-
- int GetnSurfaces();
- double Get1overP(int);
- double GetPhi(int);
- double GetTheta(int);
- double GetRadius(int);
- double GetZ(int);
- double Get1overPerror(int);
- double GetPhierror(int);
- double GetThetaerror(int);
- double GetRadiuserror(int);
- double GetZerror(int);
- double GetValue(int);
- double GetCharge();
- double LastFitMeasurement();
- double FirstFitMeasurement();
-
- const Trk::TrackParameters* GetTrackParameters(int);
-
- int GetNMeasurements();
-
- void FillVariables(const DataVector <const Trk::TrackStateOnSurface>&);
-
- void FindPotentials();
-
- GraphParameter FindParameters();
-
- GraphParameter MergeParameters(GraphParameter);
-
- private:
-
- int nmeasurements (double, double);
-
- int npotentials (double, double);
-
- double maxvalue (double, double);
-
- double coefficient_sum (std::vector<double>);
-
- double value_sum (std::vector<double>,std::vector<double>);
-
- double width_sum (std::vector<double>,std::vector<double>,std::vector<double>);
-
- //defining the tools used in this class
- ToolHandle<IMultiComponentStateCombiner> m_stateCombiner;
-
- int m_surfaceCounter;
- int m_type;
- bool m_Z_mode;
-
- mutable std::vector<double> *m_potentialvalue;
- mutable std::vector<double> *m_potentialwidth;
-
- double m_charge;
- std::vector<double> *m_1overPvalue;
- std::vector<double> *m_1overPvalueerror;
- std::vector<double> *m_PhiValue;
- std::vector<double> *m_PhiValueerror;
- std::vector<double> *m_ThetaValue;
- std::vector<double> *m_ThetaValueerror;
- std::vector<double> *m_RadiusValue;
- std::vector<double> *m_Zvalue;
- std::vector<double> *m_FitValue;
- std::vector<const Trk::TrackParameters*> *m_trackParameters;
-
- GraphParameter m_graphparameter;
- GraphParameter m_mergedparameter;
-
- };
-
-}
-#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
index faae31bfe1fa7ace7af3455acd4ba8309f373bfc..4b91eb5328b6f7e3808198fb996358344c1e995e 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/components/TrkGaussianSumFilter_entries.cxx
@@ -20,8 +20,6 @@
#include "TrkGaussianSumFilter/GsfOutlierLogic.h"
#include "TrkGaussianSumFilter/GsfExtrapolator.h"
#include "TrkGaussianSumFilter/GsfSmoother.h"
-#include "TrkGaussianSumFilter/BremFind.h"
-//#include "TrkDetDescrInterfaces/ITrackingGeometrySvc.h"
DECLARE_COMPONENT( Trk::QuickCloseComponentsMultiStateMerger )
DECLARE_COMPONENT( Trk::MultiComponentStateModeCalculator )
@@ -45,5 +43,4 @@ DECLARE_COMPONENT( Trk::ForwardGsfFitter )
DECLARE_COMPONENT( Trk::GsfOutlierLogic )
DECLARE_COMPONENT( Trk::GsfExtrapolator )
DECLARE_COMPONENT( Trk::GsfSmoother )
-DECLARE_COMPONENT( Trk::BremFind )