diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
index f12ccaa7d336bce3f9be4f0758d4ca06e1b878b9..2481c6f2eeb0109ac08c55a6e133e341e6645268 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
@@ -408,6 +408,7 @@ SegmentFitAlg::checkFit(const Tracker::SegmentFitAlg::fitInfo& seed) const
return false;
}
if (m_tanThetaCut > 0.0 && pow(seed.fitParams(2),2.0) + pow(seed.fitParams(3), 2.0) > m_tanThetaCut * m_tanThetaCut) return false;
+ if (m_tanThetaXZCut > 0.0 && pow(seed.fitParams(2), 2.0) > m_tanThetaXZCut * m_tanThetaXZCut) return false;
for (auto i : seed.candidates)
{
double z = seed.clusters[i]->cluster.globalPosition().z();
@@ -518,7 +519,7 @@ SegmentFitAlg::selectFits(const FitMap& fits) const
// smaller than m_minClustersPerFit
info.remove_if([&](std::shared_ptr<fitInfo> p) {return
((p->clusterMask & selected->clusterMask).count() > m_sharedHitFraction * p->clusterMask.count())
- || (p->clusterMask.count() < m_minClustersPerFit) ;}
+ || ((selectedFits.size() >= 4) && (p->clusterMask.count() < m_minClustersPerFit)) ;}
);
}
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
index 757eb1a22752acbe72e0da7a0c6c7a89878bf76a..1b9058ed9aa7e9a1a59e3110e1b4cf136aca910d 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
@@ -413,6 +413,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
DoubleProperty m_waferTolerance { this, "WaferTolerance", 3 * Gaudi::Units::mm, "Wafer boundary tolerance for extrapolated fit." };
DoubleProperty m_yResidualCut { this, "ResidualCut", 3 * Gaudi::Units::mm, "Cluster y tolerance for compatibility with extrapolated fit." };
DoubleProperty m_tanThetaCut { this, "TanThetaCut", 0.0, "Maximum accepted tangent of track angle from beam axis; 0 means no cut."};
+ DoubleProperty m_tanThetaXZCut { this, "TanThetaXZCut", 0.0, "Maximum accepted tangent of track angle from beam axis in x direction; 0 means no cut."};
DoubleProperty m_reducedChi2Cut { this, "ReducedChi2Cut", 10.0, "Maximum accepted chi^2 per degree of freedom for final fits; 0 means no cut." };
DoubleProperty m_sharedHitFraction { this, "SharedHitFraction", -1., "Fraction of hits which are allowed to be shared between two fits." };
IntegerProperty m_minClustersPerFit { this, "MinClustersPerFit", 4, "Minimum number of clusters a fit has to have." };
diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
index c5c0c5e454ab0b43622689bda4f7c3691ee1abca..22f075abb55848488419bc7129574a6a6a6741fe 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
@@ -27,6 +27,7 @@ atlas_add_library( FaserActsKalmanFilterLib
atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/ActsTrackSeedTool.h
FaserActsKalmanFilter/CircleFit.h
+ FaserActsKalmanFilter/CircleFitTrackSeedTool.h
FaserActsKalmanFilter/CKF2.h
FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
FaserActsKalmanFilter/EffPlotTool.h
@@ -34,12 +35,14 @@ atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/FaserActsGeometryContainers.h
FaserActsKalmanFilter/FaserActsKalmanFilterAlg.h
FaserActsKalmanFilter/FaserActsRecMultiTrajectory.h
+ FaserActsKalmanFilter/GhostBusters.h
FaserActsKalmanFilter/MyTrackSeedTool.h
FaserActsKalmanFilter/IdentifierLink.h
FaserActsKalmanFilter/IndexSourceLink.h
FaserActsKalmanFilter/ITrackFinderTool.h
FaserActsKalmanFilter/ITrackSeedTool.h
FaserActsKalmanFilter/KalmanFitterTool.h
+ FaserActsKalmanFilter/LinearFit.h
# FaserActsKalmanFilter/ClusterTrackSeedTool.h
# FaserActsKalmanFilter/TruthTrackFinderTool.h
FaserActsKalmanFilter/Measurement.h
@@ -49,6 +52,7 @@ atlas_add_component(FaserActsKalmanFilter
FaserActsKalmanFilter/ResPlotTool.h
FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
+ FaserActsKalmanFilter/SeedingAlg.h
# FaserActsKalmanFilter/SegmentFitClusterTrackFinderTool.h
# FaserActsKalmanFilter/SegmentFitTrackFinderTool.h
FaserActsKalmanFilter/SimWriterTool.h
@@ -64,17 +68,21 @@ atlas_add_component(FaserActsKalmanFilter
# FaserActsKalmanFilter/TruthSeededTrackFinderTool.h
FaserActsKalmanFilter/ThreeStationTrackSeedTool.h
src/ActsTrackSeedTool.cxx
+ src/CircleFit.cxx
+ src/CircleFitTrackSeedTool.cxx
src/CKF2.cxx
# src/ClusterTrackSeedTool.cxx
src/CombinatorialKalmanFilterAlg.cxx
src/EffPlotTool.cxx
src/FaserActsKalmanFilterAlg.cxx
+ src/GhostBusters.cxx
src/MyTrackSeedTool.cxx
src/KalmanFitterTool.cxx
# src/MultiTrackFinderTool.cxx
src/PerformanceWriterTool.cxx
src/PlotHelpers.cxx
src/ResPlotTool.cxx
+ src/SeedingAlg.cxx
src/RootTrajectoryStatesWriterTool.cxx
src/RootTrajectorySummaryWriterTool.cxx
# src/SegmentFitClusterTrackFinderTool.cxx
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
index 7e3b0e76af52738963d42f7cd716a1c32644ca47..88c659754820b2309da641d6c711b937cd299201 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
@@ -131,6 +131,7 @@ private:
Gaudi::Property<double> m_maxSteps {this, "maxSteps", 10000};
Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15};
Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10};
+ Gaudi::Property<double> m_origin {this, "origin", 0};
Gaudi::Property<size_t> m_nTrajectories {this, "nTrajectories", 2};
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
index 664b7005ddfdb31998dc3a5665fdd871a9d8ce31..780a476a722707f280b427a1678f16fcfb5b1cdf 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFit.h
@@ -5,6 +5,25 @@
#include <cmath>
#include <vector>
+namespace CircleFit {
+
+
+class CircleData {
+public:
+ CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint *> &spacePoints);
+ CircleData(const std::vector<Amg::Vector3D> &spacePoints);
+ double meanX() const;
+ double meanY() const;
+ double x(int i) const;
+ double y(int i) const;
+ int size() const;
+
+private:
+ std::vector<double> m_x{};
+ std::vector<double> m_y{};
+ int m_size = 0;
+};
+
struct Circle {
public:
@@ -19,170 +38,9 @@ public:
double j = 0;
};
+double sigma(const CircleData &data, const Circle &circle);
+Circle circleFit(const CircleData &data);
-class CircleData {
-public:
- CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint*> &spacePoints) {
- for (auto sp : spacePoints) {
- m_x.push_back(sp->globalPosition().z());
- m_y.push_back(sp->globalPosition().y());
- }
- m_size = spacePoints.size();
- }
-
- double meanX() const {
- double mean = 0;
- for (double i : m_x) mean += i;
- return mean / m_size;
- }
-
- double meanY() const {
- double mean = 0;
- for (double i : m_y) mean += i;
- return mean / m_size;
- }
-
- double x(int i) const {
- return m_x[i];
- }
-
- double y(int i) const {
- return m_y[i];
- }
-
- int size() const {
- return m_size;
- }
-
-private:
- std::vector<double> m_x {};
- std::vector<double> m_y {};
- int m_size = 0;
-};
-
-
-double Sigma (CircleData& data, Circle& circle) {
- double sum=0.,dx,dy;
-
- for (int i=0; i<data.size(); i++) {
- dx = data.x(i) - circle.cx;
- dy = data.y(i) - circle.cy;
- sum += (sqrt(dx*dx+dy*dy) - circle.r) * (sqrt(dx*dx+dy*dy) - circle.r);
- }
- return sqrt(sum/data.size());
-}
-
-
-Circle CircleFit(CircleData& data)
-/*
- Circle fit to a given set of data points (in 2D)
-
- This is an algebraic fit, due to Taubin, based on the journal article
-
- G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
- Space Curves Defined By Implicit Equations, With
- Applications To Edge And Range Image Segmentation",
- IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
-
- Input: data - the class of data (contains the given points):
-
- data.size() - the number of data points
- data.X[] - the array of X-coordinates
- data.Y[] - the array of Y-coordinates
-
- Output:
- circle - parameters of the fitting circle:
-
- circle.a - the X-coordinate of the center of the fitting circle
- circle.b - the Y-coordinate of the center of the fitting circle
- circle.r - the radius of the fitting circle
- circle.s - the root mean square error (the estimate of sigma)
- circle.j - the total number of iterations
-
- The method is based on the minimization of the function
-
- sum [(x-a)^2 + (y-b)^2 - R^2]^2
- F = -------------------------------
- sum [(x-a)^2 + (y-b)^2]
-
- This method is more balanced than the simple Kasa fit.
-
- It works well whether data points are sampled along an entire circle or
- along a small arc.
-
- It still has a small bias and its statistical accuracy is slightly
- lower than that of the geometric fit (minimizing geometric distances),
- but slightly higher than that of the very similar Pratt fit.
- Besides, the Taubin fit is slightly simpler than the Pratt fit
-
- It provides a very good initial guess for a subsequent geometric fit.
-
- Nikolai Chernov (September 2012)
-
-*/
-{
- int i,iter,IterMAX=99;
-
- double Xi,Yi,Zi;
- double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
- double A0,A1,A2,A22,A3,A33;
- double Dy,xnew,x,ynew,y;
- double DET,Xcenter,Ycenter;
-
- Circle circle;
-
- Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
- for (i=0; i<data.size(); i++) {
- Xi = data.x(i) - data.meanX();
- Yi = data.y(i) - data.meanY();
- Zi = Xi*Xi + Yi*Yi;
-
- Mxy += Xi*Yi;
- Mxx += Xi*Xi;
- Myy += Yi*Yi;
- Mxz += Xi*Zi;
- Myz += Yi*Zi;
- Mzz += Zi*Zi;
- }
- Mxx /= data.size();
- Myy /= data.size();
- Mxy /= data.size();
- Mxz /= data.size();
- Myz /= data.size();
- Mzz /= data.size();
-
- Mz = Mxx + Myy;
- Cov_xy = Mxx*Myy - Mxy*Mxy;
- Var_z = Mzz - Mz*Mz;
- A3 = 4*Mz;
- A2 = -3*Mz*Mz - Mzz;
- A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
- A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
- A22 = A2 + A2;
- A33 = A3 + A3 + A3;
-
- for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) {
- Dy = A1 + x*(A22 + A33*x);
- xnew = x - y/Dy;
- if ((xnew == x)||(!std::isfinite(xnew))) break;
- ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
- if (abs(ynew)>=abs(y)) break;
- x = xnew; y = ynew;
- }
-
- DET = x*x - x*Mz + Cov_xy;
- Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
- Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
-
- circle.cx = Xcenter + data.meanX();
- circle.cy = Ycenter + data.meanY();
- circle.r = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
- circle.s = Sigma(data,circle);
- circle.i = 0;
- circle.j = iter;
-
- return circle;
-}
-
+} // namespace CircleFit
-#endif // FASERACTSKALMANFILTER_CIRCLEFIT_H
\ No newline at end of file
+#endif // FASERACTSKALMANFILTER_CIRCLEFIT_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..51fcd477af77c1f576db1954921d36dc62758cba
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CircleFitTrackSeedTool.h
@@ -0,0 +1,160 @@
+#ifndef FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
+#define FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
+
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "Gaudi/Property.h"
+#include "GaudiKernel/IInterface.h"
+#include "GaudiKernel/StatusCode.h"
+
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsKalmanFilter/ITrackSeedTool.h"
+#include "TrkTrack/TrackCollection.h"
+// #include "TrackerSimData/TrackerSimDataCollection.h"
+// #include "GeneratorObjects/McEventCollection.h"
+#include <array>
+#include <memory>
+#include <string>
+#include <vector>
+#include <boost/dynamic_bitset.hpp>
+
+typedef boost::dynamic_bitset<> ClusterSet;
+
+class FaserSCT_ID;
+namespace TrackerDD { class SCT_DetectorManager; }
+
+class CircleFitTrackSeedTool : public extends<AthAlgTool, ITrackSeedTool> {
+public:
+ CircleFitTrackSeedTool(const std::string& type, const std::string& name, const IInterface* parent);
+ virtual ~CircleFitTrackSeedTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+ virtual StatusCode run() override;
+
+ const std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> initialTrackParameters() const override;
+ const std::shared_ptr<const Acts::Surface> initialSurface() const override;
+ const std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks() const override;
+ const std::shared_ptr<IdentifierLink> idLinks() const override;
+ const std::shared_ptr<std::vector<Measurement>> measurements() const override;
+ const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters() const override;
+ const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters() const override;
+ const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints() const override;
+
+private:
+ struct Segment {
+ public:
+ Segment(const Trk::Track* track, const FaserSCT_ID *idHelper);
+ int station;
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters;
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints;
+ ClusterSet clusterSet;
+ Acts::Vector3 position;
+ std::vector<Acts::Vector3> fakePositions;
+ Acts::Vector3 momentum;
+ };
+
+ struct Seed {
+ Seed(const std::vector<Segment> &segments);
+
+ ClusterSet clusterSet;
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters;
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints;
+ std::vector<Acts::Vector3> positions;
+ std::vector<Acts::Vector3> fakePositions;
+
+ double c0, c1, cx, cy, r, chi2, momentum, charge;
+ size_t size;
+ Acts::CurvilinearTrackParameters get_params(double origin, Acts::BoundSymMatrix cov) const;
+
+ private:
+ void getChi2();
+ double getX(double z);
+ double getY(double z);
+ void fit();
+ void fakeFit(double B=0.55);
+ void linearFit();
+ double dx, dy;
+ double m_MeV2GeV = 1e-3;
+ double m_sigma_x = 0.8;
+ double m_sigma_y = 0.016;
+ };
+
+ void go(const std::array<std::vector<Segment>, 4> &v, std::vector<Segment> &combination,
+ std::vector<Seed> &seeds, int offset, int k);
+
+ static std::map<Identifier, Index> indexMap;
+ static std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> spacePointMap;
+
+ std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> m_initialTrackParameters;
+ std::shared_ptr<const Acts::Surface> m_initialSurface;
+ std::shared_ptr<std::vector<IndexSourceLink>> m_sourceLinks {};
+ std::shared_ptr<IdentifierLink> m_idLinks {};
+ std::shared_ptr<std::vector<Measurement>> m_measurements {};
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> m_clusters {};
+ std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> m_seedClusters {};
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> m_spacePoints {};
+
+ const FaserSCT_ID* m_idHelper {nullptr};
+ const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
+
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool { this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "SegmentFit", "Input track collection name" };
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainerKey { this, "ClusterContainer", "SCT_ClusterContainer"};
+ SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey { this, "SpacePoints", "SCT_SpacePointContainer"};
+ // SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey { this, "McEventCollection", "TruthEvent"};
+ // SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey { this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+
+ Gaudi::Property<double> m_std_cluster {this, "std_cluster", 0.04};
+ // covariance of the initial parameters
+ Gaudi::Property<double> m_covLoc0 {this, "covLoc0", 1};
+ Gaudi::Property<double> m_covLoc1 {this, "covLoc1", 1};
+ Gaudi::Property<double> m_covPhi {this, "covPhi", 1};
+ Gaudi::Property<double> m_covTheta {this, "covTheta", 1};
+ Gaudi::Property<double> m_covQOverP {this, "covQOverP", 1};
+ Gaudi::Property<double> m_covTime {this, "covTime", 1};
+
+ Gaudi::Property<double> m_origin {this, "origin", 0, "z position of the reference surface"};
+};
+
+inline const std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>>
+CircleFitTrackSeedTool::initialTrackParameters() const {
+ return m_initialTrackParameters;
+}
+
+inline const std::shared_ptr<const Acts::Surface>
+CircleFitTrackSeedTool::initialSurface() const {
+ return m_initialSurface;
+}
+
+inline const std::shared_ptr<std::vector<IndexSourceLink>>
+CircleFitTrackSeedTool::sourceLinks() const {
+ return m_sourceLinks;
+}
+
+inline const std::shared_ptr<IdentifierLink>
+CircleFitTrackSeedTool::idLinks() const {
+ return m_idLinks;
+}
+
+inline const std::shared_ptr<std::vector<Measurement>>
+CircleFitTrackSeedTool::measurements() const {
+ return m_measurements;
+}
+
+inline const std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>>
+CircleFitTrackSeedTool::clusters() const {
+ return m_clusters;
+}
+
+inline const std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>>
+CircleFitTrackSeedTool::seedClusters() const {
+ return m_seedClusters;
+}
+
+inline const std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>>
+CircleFitTrackSeedTool::spacePoints() const {
+ return m_spacePoints;
+}
+
+#endif // FASERACTSKALMANFILTER_CIRCLEFITTRACKSEEDTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
new file mode 100644
index 0000000000000000000000000000000000000000..7fd49edd52e43f1d5021d0ee80744bae5daa41bc
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
@@ -0,0 +1,74 @@
+#ifndef FASERACTSKALMANFILTER_GHOSTBUSTERS_H
+#define FASERACTSKALMANFILTER_GHOSTBUSTERS_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
+#include "FaserActsKalmanFilter/TrackClassification.h"
+
+
+class TTree;
+class FaserSCT_ID;
+
+class GhostBusters : public AthReentrantAlgorithm, AthHistogramming {
+public:
+ GhostBusters(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~GhostBusters() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+ struct Segment {
+ public:
+ Segment(const Trk::Track *track, const TrackerSimDataCollection &simDataCollection,
+ const FaserSCT_ID *idHelper);
+ inline int station() const { return m_station; }
+ inline Amg::Vector3D position() const { return m_position; }
+ inline double x() const { return m_position.x(); }
+ inline double y() const { return m_position.y(); }
+ inline double z() const { return m_position.z(); }
+ inline double chi2() const { return m_chi2; }
+ inline size_t size() const { return m_size; }
+ inline bool isGhost() const { return m_isGhost; }
+ inline size_t majorityHits() const { return m_particleHitCounts.empty() ? -1 : m_particleHitCounts.front().hitCount; }
+ inline const Trk::Track *track() const { return m_track; }
+ inline void setGhost() { m_isGhost = true; }
+ private:
+ std::vector<ParticleHitCount> m_particleHitCounts {};
+ size_t m_size;
+ Amg::Vector3D m_position;
+ int m_station;
+ double m_chi2;
+ bool m_isGhost = false;
+ const Trk::Track *m_track;
+ };
+
+ ServiceHandle <ITHistSvc> m_histSvc;
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputCollection", "Segments", "Output track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "SegmentFit", "Input track collection name" };
+ SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey { this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+ DoubleProperty m_xTolerance {this, "xTolerance", 0.5};
+ DoubleProperty m_yTolerance {this, "yTolerance", 0.25};
+
+ const FaserSCT_ID *m_idHelper;
+
+ mutable TTree* m_tree;
+ mutable unsigned int m_event_number;
+ mutable double m_x;
+ mutable double m_y;
+ mutable double m_z;
+ mutable double m_chi2;
+ mutable unsigned int m_station;
+ mutable unsigned int m_majorityHits;
+ mutable unsigned int m_size;
+ mutable bool m_isGhost;
+};
+
+inline const ServiceHandle <ITHistSvc> &GhostBusters::histSvc() const {
+ return m_histSvc;
+}
+
+#endif // FASERACTSKALMANFILTER_GHOSTBUSTERS_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
index 9c26faeb3f354618133ae478fca9a4f3cea947bd..1a5ba409bb30e1c62e7f79083f6599663149a828 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
@@ -59,6 +59,7 @@ private:
Gaudi::Property<bool> m_summaryWriter {this, "SummaryWriter", false};
Gaudi::Property<bool> m_statesWriter {this, "StatesWriter", false};
Gaudi::Property<bool> m_noDiagnostics {this, "noDiagnostics", true, "Set ACTS logging level to INFO and do not run performance writer, states writer or summary writer"};
+ Gaudi::Property<double> m_origin {this, "origin", 0.};
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h
new file mode 100644
index 0000000000000000000000000000000000000000..4531156216fceaab578d37662e514e9005226b5f
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/LinearFit.h
@@ -0,0 +1,30 @@
+#ifndef FASERACTSKALMANFILTER_LINEARFIT_H
+#define FASERACTSKALMANFILTER_LINEARFIT_H
+
+#include "eigen3/Eigen/Dense"
+
+namespace LinearFit {
+
+using Vector2 = Eigen::Matrix<double, 2, 1>;
+
+std::pair <Vector2, Vector2 > linearFit(const std::vector<Amg::Vector3D> &spacePoints) {
+ std::vector<Vector2> points {};
+ for (auto sp: spacePoints) {
+ points.push_back({sp.z(), sp.x()});
+ }
+ size_t nPoints = points.size();
+ Eigen::Matrix< Vector2::Scalar, Eigen::Dynamic, Eigen::Dynamic > centers(nPoints, 2);
+ for (size_t i = 0; i < nPoints; ++i) centers.row(i) = points[i];
+
+ Vector2 origin = centers.colwise().mean();
+ Eigen::MatrixXd centered = centers.rowwise() - origin.transpose();
+ Eigen::MatrixXd cov = centered.adjoint() * centered;
+ Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
+ Vector2 axis = eig.eigenvectors().col(1).normalized();
+
+ return std::make_pair(origin, axis);
+}
+
+} // namespace LinearFit
+
+#endif // FASERACTSKALMANFILTER_LINEARFIT_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..49fd1c901ab2f0159cd033401ddde589b61262e6
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SeedingAlg.h
@@ -0,0 +1,23 @@
+#ifndef FASERACTSKALMANFILTER_SEEDINGALG_H
+#define FASERACTSKALMANFILTER_SEEDINGALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "FaserActsKalmanFilter/ITrackSeedTool.h"
+#include <boost/dynamic_bitset.hpp>
+
+
+class SeedingAlg : public AthAlgorithm {
+public:
+ SeedingAlg(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~SeedingAlg() = default;
+
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+private:
+ ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "CircleFitTrackSeedTool"};
+};
+
+#endif // FASERACTSKALMANFILTER_SEEDINGALG_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
index 790cdb8f95643d270fe12e6e7a38ba0b6739382e..2ca594c7dfd08c322c08ae193239fa8cafbc5ae9 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/TrackClassification.h
@@ -15,4 +15,9 @@ void identifyContributingParticles(
const FaserActsRecMultiTrajectory& trajectories, size_t tip,
std::vector<ParticleHitCount>& particleHitCounts);
+void identifyContributingParticles(
+ const TrackerSimDataCollection& simDataCollection,
+ const std::vector<const Tracker::FaserSCT_Cluster*> clusters,
+ std::vector<ParticleHitCount>& particleHitCounts);
+
#endif // FASERACTSKALMANFILTER_TRACKCLASSIFICATION_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
index 728cc83e1cf04a37a67486bcaee6662e5dfb73ad..41550fcfe4789c7b0fbb226a43bb0b3242f3cf45 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
@@ -34,10 +34,11 @@ def CKF2Cfg(flags, **kwargs):
acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
acc.merge(acts_tracking_geometry_svc )
-# track_seed_tool = CompFactory.ClusterTrackSeedTool()
- track_seed_tool = CompFactory.ThreeStationTrackSeedTool()
+ # track_seed_tool = CompFactory.ClusterTrackSeedTool()
# track_seed_tool = CompFactory.ActsTrackSeedTool()
# track_seed_tool = CompFactory.MyTrackSeedTool()
+ # track_seed_tool = CompFactory.ThreeStationTrackSeedTool()
+ track_seed_tool = CompFactory.CircleFitTrackSeedTool()
sigma_loc0 = 1.9e-2
sigma_loc1 = 9e-1
sigma_phi = 3.3e-2
@@ -54,6 +55,7 @@ def CKF2Cfg(flags, **kwargs):
track_seed_tool.covQOverP = initial_variance_inflation[4] * sigma_qop * sigma_qop
track_seed_tool.std_cluster = 0.0231
track_seed_tool.origin = 0
+ track_seed_tool.TrackCollection = "Segments"
trajectory_states_writer_tool = CompFactory.RootTrajectoryStatesWriterTool()
trajectory_states_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
@@ -107,7 +109,7 @@ def CKF2Cfg(flags, **kwargs):
ckf.ActsLogging = "INFO"
ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool
ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool
- ckf.PerformanceWriter = trajectory_performance_writer_tool
+ ckf.PerformanceWriterTool = trajectory_performance_writer_tool
ckf.nMax = 10
ckf.chi2Max = 25
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..cd07ccd71e798da0f7668f3f4c8512956f17683d
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
@@ -0,0 +1,35 @@
+"""
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+
+def GhostBusters_OutputCfg(flags):
+ acc = ComponentAccumulator()
+ itemList = ["xAOD::EventInfo#*",
+ "xAOD::EventAuxInfo#*",
+ "TrackCollection#*",
+ ]
+ acc.merge(OutputStreamCfg(flags, "ESD", itemList))
+ ostream = acc.getEventAlgo("OutputStreamESD")
+ ostream.TakeItemsFromInput = True
+ return acc
+
+
+def GhostBustersCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+ acc.addEventAlgo(CompFactory.GhostBusters(**kwargs))
+ acc.merge(GhostBusters_OutputCfg(flags))
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST2 DATAFILE='GhostBusters.root' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..e9b200ccb1e8c209bb1b27cb1555e8b8b27ec7bc
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/SeedingConfig.py
@@ -0,0 +1,16 @@
+"""
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+
+
+def SeedingCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+ seedingAlg = CompFactory.SeedingAlg(**kwargs)
+ acc.addEventAlgo(seedingAlg)
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 355e4a930024de1239d222e24fddf1f7560b6feb..14d6b22aa987975dddc87875cc0dbb2fef7b7be2 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -98,11 +98,7 @@ StatusCode CKF2::execute() {
m_trackSeedTool->initialTrackParameters();
std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks =
m_trackSeedTool->sourceLinks();
- std::shared_ptr<IdentifierLink> idLinks = m_trackSeedTool->idLinks();
- std::vector<Identifier> ids {};
- for (const auto& idLink : *idLinks) {
- ids.push_back(idLink.second);
- }
+
std::shared_ptr<std::vector<Measurement>> measurements = m_trackSeedTool->measurements();
std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters = m_trackSeedTool->clusters();
std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints = m_trackSeedTool->spacePoints();
@@ -124,7 +120,7 @@ StatusCode CKF2::execute() {
rotation.col(0) = Acts::Vector3(0, 0, -1);
rotation.col(1) = Acts::Vector3(0, 1, 0);
rotation.col(2) = Acts::Vector3(1, 0, 0);
- Acts::Translation3 trans(0., 0., 0.);
+ Acts::Translation3 trans(0., 0., m_origin);
Acts::Transform3 trafo(rotation * trans);
initialSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(trafo);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9f1fab11cf10edac1a497c7af498e737b5bd84c6
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFit.cxx
@@ -0,0 +1,169 @@
+#include "FaserActsKalmanFilter/CircleFit.h"
+
+namespace CircleFit {
+
+CircleData::CircleData(const std::vector<const Tracker::FaserSCT_SpacePoint *> &spacePoints) {
+ for (auto sp: spacePoints) {
+ m_x.push_back(sp->globalPosition().z());
+ m_y.push_back(sp->globalPosition().y());
+ }
+ m_size = spacePoints.size();
+}
+
+CircleData::CircleData(const std::vector<Amg::Vector3D> &spacePoints) {
+ for (auto sp: spacePoints) {
+ m_x.push_back(sp.z());
+ m_y.push_back(sp.y());
+ }
+ m_size = spacePoints.size();
+}
+
+double CircleData::meanX() const {
+ double mean = 0;
+ for (double i: m_x) mean += i;
+ return mean / m_size;
+}
+
+double CircleData::meanY() const {
+ double mean = 0;
+ for (double i: m_y) mean += i;
+ return mean / m_size;
+}
+
+double CircleData::x(int i) const {
+ return m_x[i];
+}
+
+double CircleData::y(int i) const {
+ return m_y[i];
+}
+
+int CircleData::size() const {
+ return m_size;
+}
+
+
+double sigma(const CircleData &data, const Circle &circle) {
+ double sum=0.,dx,dy;
+ for (int i=0; i<data.size(); i++) {
+ dx = data.x(i) - circle.cx;
+ dy = data.y(i) - circle.cy;
+ sum += (sqrt(dx*dx+dy*dy) - circle.r) * (sqrt(dx*dx+dy*dy) - circle.r);
+ }
+ return sqrt(sum/data.size());
+}
+
+
+Circle circleFit(const CircleData &data)
+/*
+ Circle fit to a given set of data points (in 2D)
+
+ This is an algebraic fit, due to Taubin, based on the journal article
+
+ G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
+ Space Curves Defined By Implicit Equations, With
+ Applications To Edge And Range Image Segmentation",
+ IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
+
+ Input: data - the class of data (contains the given points):
+
+ data.size() - the number of data points
+ data.X[] - the array of X-coordinates
+ data.Y[] - the array of Y-coordinates
+
+ Output:
+ circle - parameters of the fitting circle:
+
+ circle.a - the X-coordinate of the center of the fitting circle
+ circle.b - the Y-coordinate of the center of the fitting circle
+ circle.r - the radius of the fitting circle
+ circle.s - the root mean square error (the estimate of sigma)
+ circle.j - the total number of iterations
+
+ The method is based on the minimization of the function
+
+ sum [(x-a)^2 + (y-b)^2 - R^2]^2
+ F = -------------------------------
+ sum [(x-a)^2 + (y-b)^2]
+
+ This method is more balanced than the simple Kasa fit.
+
+ It works well whether data points are sampled along an entire circle or
+ along a small arc.
+
+ It still has a small bias and its statistical accuracy is slightly
+ lower than that of the geometric fit (minimizing geometric distances),
+ but slightly higher than that of the very similar Pratt fit.
+ Besides, the Taubin fit is slightly simpler than the Pratt fit
+
+ It provides a very good initial guess for a subsequent geometric fit.
+
+ Nikolai Chernov (September 2012)
+
+*/
+{
+ int i, iter, IterMAX = 999;
+
+ double Xi, Yi, Zi;
+ double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
+ double A0, A1, A2, A22, A3, A33;
+ double Dy, xnew, x, ynew, y;
+ double DET, Xcenter, Ycenter;
+
+ Circle circle;
+
+ Mxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;
+ for (i = 0; i < data.size(); i++) {
+ Xi = data.x(i) - data.meanX();
+ Yi = data.y(i) - data.meanY();
+ Zi = Xi * Xi + Yi * Yi;
+
+ Mxy += Xi * Yi;
+ Mxx += Xi * Xi;
+ Myy += Yi * Yi;
+ Mxz += Xi * Zi;
+ Myz += Yi * Zi;
+ Mzz += Zi * Zi;
+ }
+ Mxx /= data.size();
+ Myy /= data.size();
+ Mxy /= data.size();
+ Mxz /= data.size();
+ Myz /= data.size();
+ Mzz /= data.size();
+
+ Mz = Mxx + Myy;
+ Cov_xy = Mxx * Myy - Mxy * Mxy;
+ Var_z = Mzz - Mz * Mz;
+ A3 = 4 * Mz;
+ A2 = -3 * Mz * Mz - Mzz;
+ A1 = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
+ A0 = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
+ A22 = A2 + A2;
+ A33 = A3 + A3 + A3;
+
+ for (x = 0., y = A0, iter = 0; iter < IterMAX; iter++) {
+ Dy = A1 + x * (A22 + A33 * x);
+ xnew = x - y / Dy;
+ if ((xnew == x) || (!std::isfinite(xnew))) break;
+ ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
+ if (abs(ynew) >= abs(y)) break;
+ x = xnew;
+ y = ynew;
+ }
+
+ DET = x * x - x * Mz + Cov_xy;
+ Xcenter = (Mxz * (Myy - x) - Myz * Mxy) / DET / 2;
+ Ycenter = (Myz * (Mxx - x) - Mxz * Mxy) / DET / 2;
+
+ circle.cx = Xcenter + data.meanX();
+ circle.cy = Ycenter + data.meanY();
+ circle.r = sqrt(Xcenter * Xcenter + Ycenter * Ycenter + Mz);
+ circle.s = sigma(data, circle);
+ circle.i = 0;
+ circle.j = iter;
+
+ return circle;
+}
+
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a9cbb92919799f2cabf183bfcba4ffe3f611d917
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -0,0 +1,307 @@
+#include "FaserActsKalmanFilter/CircleFitTrackSeedTool.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "FaserActsKalmanFilter/CircleFit.h"
+#include "FaserActsKalmanFilter/LinearFit.h"
+#include "FaserActsKalmanFilter/TrackClassification.h"
+#include <array>
+
+
+std::map<Identifier,Index> CircleFitTrackSeedTool::indexMap {};
+std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> CircleFitTrackSeedTool::spacePointMap {};
+
+CircleFitTrackSeedTool::CircleFitTrackSeedTool(
+ const std::string& type, const std::string& name, const IInterface* parent)
+ : base_class(type, name, parent) {}
+
+
+StatusCode CircleFitTrackSeedTool::initialize() {
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_detManager, "SCT"));
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_clusterContainerKey.initialize());
+ ATH_CHECK(m_spacePointContainerKey.initialize());
+ // ATH_CHECK(m_simDataCollectionKey.initialize());
+ // ATH_CHECK(m_mcEventCollectionKey.initialize());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode CircleFitTrackSeedTool::run() {
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer {m_clusterContainerKey};
+ ATH_CHECK(clusterContainer.isValid());
+
+ SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey};
+ ATH_CHECK(spacePointContainer.isValid());
+
+ // SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey};
+ // ATH_CHECK(simData.isValid());
+
+ // SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey};
+ // ATH_CHECK(mcEvents.isValid());
+
+ using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
+ std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
+ // std::map<int, const HepMC::GenParticle*> particles {};
+ // for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ // particles[particle->barcode()] = particle;
+ // }
+
+ // std::map<Identifier, const Tracker::FaserSCT_SpacePoint*> spacePointMap {};
+ std::vector<const Tracker::FaserSCT_SpacePoint*> spacePoints {};
+ for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+ for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
+ spacePoints.push_back(spacePoint);
+ Identifier id1 = spacePoint->cluster1()->identify();
+ CircleFitTrackSeedTool::spacePointMap[id1] = spacePoint;
+ Identifier id2 = spacePoint->cluster2()->identify();
+ CircleFitTrackSeedTool::spacePointMap[id2] = spacePoint;
+ }
+ }
+
+
+ using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, 1>;
+ std::array<Acts::BoundIndices, 1> Indices = {Acts::eBoundLoc0};
+ std::vector<IndexSourceLink> sourceLinks {};
+ std::vector<Measurement> measurements {};
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters {};
+ for (const Tracker::FaserSCT_ClusterCollection* clusterCollection : *clusterContainer) {
+ for (const Tracker::FaserSCT_Cluster* cluster : *clusterCollection) {
+ clusters.push_back(cluster);
+ Identifier id = cluster->detectorElement()->identify();
+ CircleFitTrackSeedTool::indexMap[cluster->identify()] = measurements.size();
+ if (identifierMap->count(id) != 0) {
+ Acts::GeometryIdentifier geoId = identifierMap->at(id);
+ IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
+ // create measurement
+ const auto& par = cluster->localPosition();
+ Eigen::Matrix<double, 1, 1> pos {par.x(),};
+ Eigen::Matrix<double, 1, 1> cov {0.0231 * 0.0231,};
+ ThisMeasurement meas(sourceLink, Indices, pos, cov);
+ sourceLinks.push_back(sourceLink);
+ measurements.emplace_back(std::move(meas));
+ }
+ }
+ }
+
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection};
+ ATH_CHECK(trackCollection.isValid());
+
+ std::array<std::vector<Segment>, 4> segments {};
+ for (const Trk::Track* track : *trackCollection) {
+ auto s = Segment(track, m_idHelper);
+ segments[s.station].push_back(s);
+ }
+
+ std::vector<Segment> combination {};
+ std::vector<Seed> seeds {};
+ // create seeds from four stations
+ // go(segments, combination, seeds, 0, 4);
+ // create seeds from three stations
+ go(segments, combination, seeds, 0, 3);
+ // create seeds from two stations
+ if (seeds.size() < 2) {
+ go(segments, combination, seeds, 0, 2);
+ }
+
+ std::list<Seed> allSeeds;
+ for (const Seed &seed : seeds) allSeeds.push_back(seed);
+
+ allSeeds.sort([](const Seed &left, const Seed &right) {
+ if (left.size > right.size) return true;
+ if (left.size < right.size) return false;
+ if (left.chi2 < right.chi2) return true;
+ else return false;
+ });
+
+ std::vector<Seed> selectedSeeds {};
+ while (not allSeeds.empty()) {
+ Seed selected = allSeeds.front();
+ selectedSeeds.push_back(selected);
+ allSeeds.remove_if([&](const Seed &p) {
+ return (p.size <= 12) || ((p.clusterSet & selected.clusterSet).count() > 6);
+ });
+ }
+
+
+ Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
+ cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = m_covLoc0;
+ cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = m_covLoc1;
+ cov(Acts::eBoundPhi, Acts::eBoundPhi) = m_covPhi;
+ cov(Acts::eBoundTheta, Acts::eBoundTheta) = m_covTheta;
+ cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = m_covQOverP;
+ cov(Acts::eBoundTime, Acts::eBoundTime) = m_covTime;
+
+ std::vector<Acts::CurvilinearTrackParameters> initParams {};
+ for (const Seed &seed : selectedSeeds) {
+ initParams.push_back(seed.get_params(m_origin, cov));
+ }
+
+ m_initialTrackParameters = std::make_shared<std::vector<Acts::CurvilinearTrackParameters>>(initParams);
+ m_sourceLinks = std::make_shared<std::vector<IndexSourceLink>>(sourceLinks);
+ // m_idLinks = std::make_shared<IdentifierLink>(identifierLinkMap);
+ m_measurements = std::make_shared<std::vector<Measurement>>(measurements);
+ m_initialSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3 {0, 0, m_origin}, Acts::Vector3{0, 0, -1});
+ m_clusters = std::make_shared<std::vector<const Tracker::FaserSCT_Cluster*>>(clusters);
+ m_spacePoints = std::make_shared<std::vector<const Tracker::FaserSCT_SpacePoint*>>(spacePoints);
+
+ /*
+ std::cout.precision(17);
+ for (auto &seed : selectedSeeds) {
+ std::cout << "np.array([";
+ for (const Acts::Vector3 &pos : seed.fakePositions) {
+ std::cout << "[" << pos.x() << ", " << pos.y() << ", " << pos.z() << "], ";
+ }
+ std::cout << "])" << std::endl;
+ std::cout << "chi2: " << seed.chi2 << ", momentum: " << seed.momentum << ", charge: " << seed.charge << std::endl;
+ std::cout << "fit = np.array([" << seed.c1 << ", " << seed.c0 << ", " << seed.cy << ", " << seed.cx << ", " << seed.r << "])" << std::endl;
+
+ std::vector<ParticleHitCount> particleHitCounts;
+ identifyContributingParticles(*simData, seed.clusters, particleHitCounts);
+ auto ip = particles.find(particleHitCounts.front().particleId);
+ if (ip != particles.end()) {
+ const HepMC::GenParticle* particle = ip->second;
+ HepMC::FourVector momentum = particle->momentum();
+ std::cout << "true momentum: " << momentum.rho() * 0.001 << std::endl;
+ }
+ for (const ParticleHitCount &hitCount : particleHitCounts) {
+ std::cout << hitCount.particleId << " : " << hitCount.hitCount << std::endl;
+ }
+ }
+ */
+
+ indexMap.clear();
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode CircleFitTrackSeedTool::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+
+void CircleFitTrackSeedTool::go(const std::array<std::vector<Segment>, 4> &v,
+ std::vector<Segment> &combination,
+ std::vector<Seed> &seeds,
+ int offset, int k) {
+ if (k == 0) {
+ seeds.push_back(Seed(combination));
+ return;
+ }
+ for (int i = offset; i < v.size() + 1 - k; ++i) {
+ for (const auto& ve : v[i]) {
+ combination.push_back(ve);
+ go(v, combination, seeds, i+1, k-1);
+ combination.pop_back();
+ }
+ }
+}
+
+CircleFitTrackSeedTool::Segment::Segment(const Trk::Track* track, const FaserSCT_ID *idHelper) :
+ clusterSet(CircleFitTrackSeedTool::indexMap.size()) {
+ for (const Trk::TrackStateOnSurface* trackState : *(track->trackStateOnSurfaces())) {
+ auto clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*> (trackState->measurementOnTrack());
+ if (clusterOnTrack) {
+ const Tracker::FaserSCT_Cluster* cluster = clusterOnTrack->prepRawData();
+ Identifier id = cluster->identify();
+ clusters.push_back(cluster);
+ if (CircleFitTrackSeedTool::spacePointMap.count(id) > 0) {
+ const Tracker::FaserSCT_SpacePoint *sp = CircleFitTrackSeedTool::spacePointMap.at(cluster->identify());
+ if (std::find(spacePoints.begin(), spacePoints.end(), sp) == spacePoints.end()) {
+ spacePoints.push_back(sp);
+ }
+ }
+ station = idHelper->station(id);
+ clusterSet.set(CircleFitTrackSeedTool::indexMap.at(id));
+ auto fitParameters = track->trackParameters()->front();
+ position = fitParameters->position();
+ momentum = fitParameters->momentum();
+ }
+ }
+ fakePositions.push_back(position);
+ fakePositions.push_back(position - 30 * momentum.normalized());
+ fakePositions.push_back(position + 30 * momentum.normalized());
+}
+
+CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
+ clusterSet(CircleFitTrackSeedTool::indexMap.size()) {
+ for (const Segment &seg : segments) {
+ clusters.insert(clusters.end(), seg.clusters.begin(), seg.clusters.end());
+ spacePoints.insert(spacePoints.end(), seg.spacePoints.begin(), seg.spacePoints.end());
+ positions.push_back(seg.position);
+ fakePositions.insert(fakePositions.end(), seg.fakePositions.begin(), seg.fakePositions.end());
+ for (size_t i = 0; i < seg.clusterSet.size(); ++i) {
+ if (seg.clusterSet[i]) clusterSet.set(i);
+ }
+ }
+ getChi2();
+ size = clusters.size();
+}
+
+
+void CircleFitTrackSeedTool::Seed::fit() {
+ CircleFit::CircleData circleData(positions);
+ CircleFit::Circle circle = CircleFit::circleFit(circleData);
+ momentum = circle.r * 0.001 * 0.3 * 0.55;
+ charge = circle.cy > 0 ? 1 : -1;
+ // std::cout << "circle: \nnp.array([" << circle.cx << ", " << circle.cy << ", " << circle.r << "])" << std::endl;
+}
+
+void CircleFitTrackSeedTool::Seed::fakeFit(double B) {
+ CircleFit::CircleData circleData(fakePositions);
+ CircleFit::Circle circle = CircleFit::circleFit(circleData);
+ cx = circle.cx;
+ cy = circle.cy;
+ r = circle.r;
+ momentum = r * 0.3 * B * m_MeV2GeV;
+ charge = circle.cy > 0 ? 1 : -1;
+}
+
+void CircleFitTrackSeedTool::Seed::linearFit() {
+ auto [origin, dir] = LinearFit::linearFit(positions);
+ c1 = dir[1]/dir[0];
+ c0 = origin[1] - origin[0] * c1;
+}
+
+
+double CircleFitTrackSeedTool::Seed::getY(double z) {
+ double sqt = std::sqrt(-cx*cx + 2*cx*z + r*r - z*z);
+ return abs(cy - sqt) < abs(cy + sqt) ? cy - sqt : cy + sqt;
+}
+
+
+double CircleFitTrackSeedTool::Seed::getX(double z) {
+ return c0 + z * c1;
+}
+
+void CircleFitTrackSeedTool::Seed::getChi2() {
+ linearFit();
+ fakeFit();
+
+ chi2 = 0;
+ for (const Acts::Vector3 &pos : fakePositions) {
+ dy = pos.y() - getY(pos.z());
+ chi2 += (dy * dy) / (m_sigma_y * m_sigma_y);
+ }
+
+ for (const Acts::Vector3 &pos : positions) {
+ dx = pos.x() - getX(pos.x());
+ chi2 += (dx * dx) / (m_sigma_x * m_sigma_x);
+ }
+}
+
+Acts::CurvilinearTrackParameters CircleFitTrackSeedTool::Seed::get_params(double origin, Acts::BoundSymMatrix cov) const {
+ Acts::Vector3 dir = positions[1] - positions[0];
+ dir = dir.normalized();
+ Acts::Vector3 pos = positions[0] - (positions[0].z() - origin)/dir.z() * dir;
+ Acts::Vector4 pos4 {pos.x(), pos.y(), pos.z(), 0};
+ return Acts::CurvilinearTrackParameters(pos4, dir, momentum, charge, cov);
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4525f42cd575b0c3eee67893c09fca0dadeea3d6
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
@@ -0,0 +1,115 @@
+#include "FaserActsKalmanFilter/GhostBusters.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+
+
+GhostBusters::GhostBusters(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator), AthHistogramming(name),
+ m_histSvc("THistSvc/THistSvc", name) {}
+
+
+StatusCode GhostBusters::initialize() {
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_simDataCollectionKey.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+
+ m_tree = new TTree("tree", "tree");
+ m_tree->Branch("event_number", &m_event_number, "event_number/I");
+ m_tree->Branch("station", &m_station, "station/I");
+ m_tree->Branch("x", &m_x, "x/D");
+ m_tree->Branch("y", &m_y, "y/D");
+ m_tree->Branch("z", &m_z, "z/D");
+ m_tree->Branch("chi2", &m_chi2, "chi2/D");
+ m_tree->Branch("majorityHits", &m_majorityHits, "majorityHits/I");
+ m_tree->Branch("size", &m_size, "size/I");
+ m_tree->Branch("isGhost", &m_isGhost, "isGhost/B");
+ ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode GhostBusters::execute(const EventContext &ctx) const {
+ m_event_number = ctx.eventID().event_number();
+
+ SG::WriteHandle outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
+ ATH_CHECK(simData.isValid());
+
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ ATH_CHECK(trackCollection.isValid());
+
+ std::map<int, std::vector<Segment>> segments {};
+ for (const Trk::Track* track : *trackCollection) {
+ auto segment = Segment(track, *simData, m_idHelper);
+ segments[segment.station()].push_back(segment);
+ }
+
+ for (const auto &station : segments) {
+ double nGoodSegments = 0;
+ std::vector<Segment> stationSegments = station.second;
+ std::sort(stationSegments.begin(), stationSegments.end(), [](const Segment &lhs, const Segment &rhs) { return lhs.y() > rhs.y(); });
+ double maxX = stationSegments.front().x();
+ double maxY = stationSegments.front().y();
+ double minX = stationSegments.back().x();
+ double minY = stationSegments.back().y();
+ double deltaY = maxY - minY;
+ double meanX = 0.5 * (minX + maxX);
+ double meanY = 0.5 * (minY + maxY);
+ double deltaX = deltaY / (2 * std::tan(0.02));
+ for (Segment &segment : stationSegments) {
+ bool isGhost = (segment.y() > meanY - m_yTolerance * deltaY) && (segment.y() < meanY + m_yTolerance * deltaY) && (
+ ((segment.x() > meanX - (1 + m_xTolerance) * deltaX) &&
+ (segment.x() < meanX - (1 - m_xTolerance) * deltaX)) ||
+ ((segment.x() > meanX + (1 - m_xTolerance) * deltaX) && (segment.x() < meanX + (1 + m_xTolerance) * deltaX)));
+ if (isGhost) segment.setGhost();
+ if (not isGhost && segment.size() >= 5) nGoodSegments++;
+ }
+ for (const Segment &segment : stationSegments) {
+ m_x = segment.x();
+ m_y = segment.y();
+ m_z = segment.z();
+ m_chi2 = segment.chi2();
+ m_station = segment.station();
+ m_size = segment.size();
+ m_majorityHits = segment.majorityHits();
+ m_isGhost = segment.isGhost();
+ if (nGoodSegments >= 2 && segment.size() == 4) m_isGhost = true;
+ m_tree->Fill();
+ if (segment.isGhost()) continue;
+ if (nGoodSegments >= 2 && segment.size() == 4) continue;
+ outputTracks->push_back(new Trk::Track(*segment.track()));
+ }
+ }
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode GhostBusters::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+
+GhostBusters::Segment::Segment(const Trk::Track *track, const TrackerSimDataCollection &simDataCollection,
+ const FaserSCT_ID *idHelper) {
+ m_track = track;
+ m_position = track->trackParameters()->front()->position();
+ m_chi2 = track->fitQuality()->chiSquared();
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters {};
+ for (const Trk::TrackStateOnSurface* trackState : *(track->trackStateOnSurfaces())) {
+ auto clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*> (trackState->measurementOnTrack());
+ if (clusterOnTrack) {
+ clusters.emplace_back(clusterOnTrack->prepRawData());
+ m_station = idHelper->station(clusterOnTrack->identify());
+ }
+ }
+ m_size = clusters.size();
+ identifyContributingParticles(simDataCollection, clusters, m_particleHitCounts);
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index 686047ae7232deb09ca207f6599666f0bd022aa3..a8e0ad41de7b668fd360265dcf5ebe5bd12178e6 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -57,7 +57,7 @@ KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- Acts::Vector3 {0, 0, 0}, Acts::Vector3{0, 0, -1});
+ Acts::Vector3 {0, 0, m_origin}, Acts::Vector3{0, 0, -1});
Acts::GeometryContext gctx = m_trackingGeometryTool->getNominalGeometryContext().context();
Acts::MagneticFieldContext mfContext = getMagneticFieldContext(ctx);
@@ -208,7 +208,7 @@ KalmanFitterTool::getParametersFromTrack(const Trk::TrackParameters *inputParame
const Acts::BoundVector& inputVector) const {
using namespace Acts::UnitLiterals;
std::shared_ptr<const Acts::Surface> pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- Acts::Vector3 {0, 0, 0}, Acts::Vector3{0, 0, -1});
+ Acts::Vector3 {0, 0, m_origin}, Acts::Vector3{0, 0, -1});
auto atlasParam = inputParameters->parameters();
auto center = inputParameters->associatedSurface().center();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0d48e4a63c750462180753e9efb39f06d15e4641
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/SeedingAlg.cxx
@@ -0,0 +1,22 @@
+#include "FaserActsKalmanFilter/SeedingAlg.h"
+
+SeedingAlg::SeedingAlg(const std::string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name, pSvcLocator) {}
+
+
+StatusCode SeedingAlg::initialize() {
+ ATH_CHECK(m_trackSeedTool.retrieve());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SeedingAlg::execute() {
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+ ATH_CHECK(m_trackSeedTool->run());
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SeedingAlg::finalize() {
+ return StatusCode::SUCCESS;
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
index 785bf2be17d9fe0c612b0250c6eae67ae801626e..6b780fc368e874f4520af430e63d1cc7c4e02575 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
@@ -45,18 +45,49 @@ void identifyContributingParticles(
if (not state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
return true;
}
+ std::vector<int> barcodes {};
// register all particles that generated this hit
- std::vector<Identifier> ids = state.uncalibrated().hit()->rdoList();
- for (const Identifier &id : ids) {
+ for (const Identifier &id : state.uncalibrated().hit()->rdoList()) {
if (simDataCollection.count(id) == 0) {
return true;
}
const auto &deposits = simDataCollection.find(id)->second.getdeposits();
for (const TrackerSimData::Deposit &deposit : deposits) {
- increaseHitCount(particleHitCounts, deposit.first->barcode());
+ int barcode = deposit.first->barcode();
+ if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
+ barcodes.push_back(barcode);
+ increaseHitCount(particleHitCounts, deposit.first->barcode());
+ }
}
}
return true;
});
sortHitCount(particleHitCounts);
}
+
+/* Identify all particles that contribute to a trajectory.
+ * If a cluster consists of multiple RDOs we check for each from which particle(s) it has been created.
+ * And if multiple particles created a RDO we increase the hit count for each of them.
+ */
+void identifyContributingParticles(
+ const TrackerSimDataCollection& simDataCollection,
+ const std::vector<const Tracker::FaserSCT_Cluster*> clusters,
+ std::vector<ParticleHitCount>& particleHitCounts) {
+ particleHitCounts.clear();
+ for (const Tracker::FaserSCT_Cluster *cluster : clusters) {
+ std::vector<int> barcodes {};
+ for (const Identifier &id : cluster->rdoList()) {
+ if (simDataCollection.count(id) == 0) continue;
+ const auto &deposits = simDataCollection.find(id)->second.getdeposits();
+ for (const TrackerSimData::Deposit &deposit : deposits) {
+ int barcode = deposit.first->barcode();
+ if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
+ barcodes.push_back(barcode);
+ increaseHitCount(particleHitCounts, deposit.first->barcode());
+ }
+ }
+ }
+ }
+ sortHitCount(particleHitCounts);
+}
+
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
index 8d56370abefe0f7399095ac3b62d5bc387ad2f35..4c23cb7febdd977df01a7a4a4de997267dbfd4c6 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
@@ -25,7 +25,9 @@
#include "FaserActsKalmanFilter/CKF2.h"
#include "FaserActsKalmanFilter/KalmanFitterTool.h"
#include "FaserActsKalmanFilter/MyTrackSeedTool.h"
-
+#include "FaserActsKalmanFilter/SeedingAlg.h"
+#include "FaserActsKalmanFilter/CircleFitTrackSeedTool.h"
+#include "FaserActsKalmanFilter/GhostBusters.h"
DECLARE_COMPONENT(FaserActsKalmanFilterAlg)
DECLARE_COMPONENT(CombinatorialKalmanFilterAlg)
@@ -50,3 +52,6 @@ DECLARE_COMPONENT(ActsTrackSeedTool)
DECLARE_COMPONENT(CKF2)
DECLARE_COMPONENT(KalmanFitterTool)
DECLARE_COMPONENT(MyTrackSeedTool)
+DECLARE_COMPONENT(SeedingAlg)
+DECLARE_COMPONENT(CircleFitTrackSeedTool)
+DECLARE_COMPONENT(GhostBusters)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py b/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
index 5d41a03a33273d3003677932920b059df871b3d5..6e65c2e148d1a890b8f363c47fa984bd1764b3a7 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/CKF2.py
@@ -11,6 +11,7 @@ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
from FaserActsKalmanFilter.CKF2Config import CKF2Cfg
log.setLevel(DEBUG)
@@ -31,8 +32,10 @@ acc.merge(PoolReadCfg(ConfigFlags))
acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
-acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.51, MinClustersPerFit=5, TanThetaCut=0.25))
-acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+acc.merge(GhostBustersCfg(ConfigFlags))
+acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=False))
+acc.getEventAlgo("CKF2").OutputLevel = DEBUG
# logging.getLogger('forcomps').setLevel(VERBOSE)
# acc.foreach_component("*").OutputLevel = VERBOSE
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py b/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py
new file mode 100644
index 0000000000000000000000000000000000000000..297f76e61fd138c325b67e98d0ea195e062eb5f3
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/GhostBusters.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+
+import sys
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.SeedingConfig import SeedingCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "ghosts.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+# ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+# acc.getEventAlgo("Tracker::SegmentFitAlg").OutputLevel = VERBOSE
+acc.merge(GhostBustersCfg(ConfigFlags, xTolerance=0.5, yTolerance=0.5))
+acc.getEventAlgo("GhostBusters").OutputLevel = DEBUG
+
+# logging.getLogger('forcomps').setLevel(VERBOSE)
+# acc.foreach_component("*").OutputLevel = VERBOSE
+# acc.foreach_component("*ClassID*").OutputLevel = VERBOSE
+# acc.getService("StoreGateSvc").Dump = True
+# acc.getService("ConditionStore").Dump = True
+# acc.printConfig(withDetails=True)
+# ConfigFlags.dump()
+
+sc = acc.run(maxEvents=-1)
+sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py b/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py
new file mode 100644
index 0000000000000000000000000000000000000000..01e04b803a6f1ed60ddd54db02fcea952ad8849c
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/SeedingDbg.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+
+import sys
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.SeedingConfig import SeedingCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "circleFitSeeding.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.FaserVersion = "FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+# ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.51, MinClustersPerFit=5, TanThetaCut=0.25))
+acc.merge(SeedingCfg(ConfigFlags))
+acc.getEventAlgo("SeedingAlg").OutputLevel = VERBOSE
+
+# logging.getLogger('forcomps').setLevel(VERBOSE)
+# acc.foreach_component("*").OutputLevel = VERBOSE
+# acc.foreach_component("*ClassID*").OutputLevel = VERBOSE
+# acc.getService("StoreGateSvc").Dump = True
+# acc.getService("ConditionStore").Dump = True
+# acc.printConfig(withDetails=True)
+# ConfigFlags.dump()
+
+sc = acc.run(maxEvents=10)
+sys.exit(not sc.isSuccess())