diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
index 88c659754820b2309da641d6c711b937cd299201..7abe5e97b19d93c82304a8a2698822caffc37aa5 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CKF2.h
@@ -141,6 +141,7 @@ private:
ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
ToolHandle<KalmanFitterTool> m_kalmanFitterTool1 {this, "KalmanFitterTool1", "KalmanFitterTool"};
ToolHandle<KalmanFitterTool> m_kalmanFitterTool2 {this, "KalmanFitterTool2", "KalmanFitterTool"};
+ Gaudi::Property<bool> m_isMC {this, "isMC", false};
std::unique_ptr<Trk::Track> makeTrack(Acts::GeometryContext& tgContext, TrackFitterResult& fitResult) const;
std::unique_ptr<Trk::Track> makeTrack(const Acts::GeometryContext &geoCtx, const FaserActsRecMultiTrajectory &traj) const;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
index ba983f245fdce513858c27cb0a814b48313eac0b..28058887c9b6f05d384962206c8cca1a64a9d66c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h
@@ -114,6 +114,7 @@ public:
Gaudi::Property<double> m_maxSteps {this, "maxSteps", 1000};
Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15};
Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10};
+ Gaudi::Property<bool> m_isMC {this, "isMC", false};
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/GhostBusters.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
index 7fd49edd52e43f1d5021d0ee80744bae5daa41bc..f8ac55d6cbb62c20ac2de7ce9c8b5741a2bfe8cb 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/GhostBusters.h
@@ -23,8 +23,9 @@ public:
private:
struct Segment {
public:
- Segment(const Trk::Track *track, const TrackerSimDataCollection &simDataCollection,
- const FaserSCT_ID *idHelper);
+ Segment(const Trk::Track *track, const FaserSCT_ID *idHelper);
+ // 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(); }
@@ -33,11 +34,11 @@ private:
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 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 {};
+ // std::vector<ParticleHitCount> m_particleHitCounts {};
size_t m_size;
Amg::Vector3D m_position;
int m_station;
@@ -49,7 +50,7 @@ private:
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"};
+ // SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey { this, "TrackerSimDataCollection", "SCT_SDO_Map"};
DoubleProperty m_xTolerance {this, "xTolerance", 0.5};
DoubleProperty m_yTolerance {this, "yTolerance", 0.25};
@@ -62,7 +63,7 @@ private:
mutable double m_z;
mutable double m_chi2;
mutable unsigned int m_station;
- mutable unsigned int m_majorityHits;
+ // mutable unsigned int m_majorityHits;
mutable unsigned int m_size;
mutable bool m_isGhost;
};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
index 1a5ba409bb30e1c62e7f79083f6599663149a828..efca583869ebb4cf60bd1458ddc913abb5f0b1d3 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/KalmanFitterTool.h
@@ -42,7 +42,7 @@ public:
virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const;
std::unique_ptr<Trk::Track> fit(const EventContext& ctx, const Trk::Track& inputTrack,
std::vector<FaserActsRecMultiTrajectory> &trajectories,
- const Acts::BoundVector& inputVector) const;
+ const Acts::BoundVector& inputVector, bool isMC) const;
private:
const FaserSCT_ID* m_idHelper {nullptr};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
index ee8ef5d25a9f808d4ed8ae87787b437188969b49..6f0885d33ec1200edbd28761352a36840fd11e47 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
@@ -29,7 +29,7 @@ public:
StatusCode initialize() override;
StatusCode finalize() override;
- StatusCode write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories) const;
+ StatusCode write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const;
private:
SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey {this, "McEventCollection", "TruthEvent"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
index 33b5ee882b8a5f4358cb8b65203704f4bbd9c9a3..90cf6cf84ca8a95b41397bd00dea30d85135967f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h
@@ -25,7 +25,7 @@ class RootTrajectorySummaryWriterTool : public AthAlgTool {
public:
RootTrajectorySummaryWriterTool(const std::string& type, const std::string& name, const IInterface* parent);
~RootTrajectorySummaryWriterTool() override = default;
- StatusCode write(const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories) const;
+ StatusCode write( const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories, bool isMC) const;
StatusCode initialize() override;
StatusCode finalize() override;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..f6e894e19d33d69d693e16c6f768b8c3ec48eec8
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/TI12CKF2Config.py
@@ -0,0 +1,123 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations
+
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+
+
+def FaserActsAlignmentCondAlgCfg(flags, **kwargs):
+ acc = ComponentAccumulator()
+ acc.addCondAlgo(CompFactory.FaserActsAlignmentCondAlg(name="FaserActsAlignmentCondAlg", **kwargs))
+ return acc
+
+
+def CKF2_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 TI12CKF2Cfg(flags, **kwargs):
+ # acc = ComponentAccumulator()
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+ acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(acts_tracking_geometry_svc )
+
+ # 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
+ sigma_theta = 2.9e-4
+ p = 1000
+ sigma_p = 0.1 * p
+ sigma_qop = sigma_p / (p * p)
+ # initial_variance_inflation = [1000, 1000, 100, 100, 10000]
+ initial_variance_inflation = [100, 100, 100, 100, 1000]
+ track_seed_tool.covLoc0 = initial_variance_inflation[0] * sigma_loc1 * sigma_loc1
+ track_seed_tool.covLoc1 = initial_variance_inflation[1] * sigma_loc0 * sigma_loc0
+ track_seed_tool.covPhi = initial_variance_inflation[2] * sigma_phi * sigma_phi
+ track_seed_tool.covTheta = initial_variance_inflation[3] * sigma_theta * sigma_theta
+ track_seed_tool.covQOverP = initial_variance_inflation[4] * sigma_qop * sigma_qop
+ track_seed_tool.std_cluster = 0.0231
+ origin = -1900
+ track_seed_tool.origin = origin
+ track_seed_tool.TrackCollection = "Segments"
+
+ trajectory_states_writer_tool = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_states_writer_tool1 = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool1.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_states_writer_tool1.FilePath = "track_states_ckf1.root"
+ trajectory_states_writer_tool2 = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool2.FilePath = "track_states_ckf2.root"
+ trajectory_states_writer_tool2.noDiagnostics = kwargs["noDiagnostics"]
+
+ trajectory_summary_writer_tool = CompFactory.RootTrajectorySummaryWriterTool(**kwargs)
+ trajectory_summary_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_summary_writer_tool1 = CompFactory.RootTrajectorySummaryWriterTool()
+ trajectory_summary_writer_tool1.FilePath = "track_summary_ckf1.root"
+ trajectory_summary_writer_tool1.noDiagnostics = kwargs["noDiagnostics"]
+ trajectory_summary_writer_tool2 = CompFactory.RootTrajectorySummaryWriterTool(**kwargs)
+ trajectory_summary_writer_tool2.FilePath = "track_summary_ckf2.root"
+ trajectory_summary_writer_tool2.noDiagnostics = kwargs["noDiagnostics"]
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ # trajectory_performance_writer_tool = CompFactory.PerformanceWriterTool("PerformanceWriterTool", **kwargs)
+ # trajectory_performance_writer_tool.ExtrapolationTool = actsExtrapolationTool
+ # trajectory_performance_writer_tool.noDiagnostics = kwargs["noDiagnostics"]
+
+
+ ckf = CompFactory.CKF2(**kwargs)
+ kalman_fitter1 = CompFactory.KalmanFitterTool(name="fitterTool1", **kwargs)
+ kalman_fitter1.noDiagnostics = kwargs["noDiagnostics"]
+ kalman_fitter1.ActsLogging = "INFO"
+ kalman_fitter1.SummaryWriter = True
+ kalman_fitter1.StatesWriter = True
+ kalman_fitter1.SeedCovarianceScale = 10
+ kalman_fitter1.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool1
+ kalman_fitter1.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool1
+ kalman_fitter1.origin = origin
+ ckf.KalmanFitterTool1 = kalman_fitter1
+
+ kalman_fitter2 = CompFactory.KalmanFitterTool(name="fitterTool2", **kwargs)
+ kalman_fitter2.noDiagnostics = kwargs["noDiagnostics"]
+ kalman_fitter2.ActsLogging = "INFO"
+ kalman_fitter2.SummaryWriter = True
+ kalman_fitter2.StatesWriter = True
+ kalman_fitter2.SeedCovarianceScale = 1
+ kalman_fitter2.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool2
+ kalman_fitter2.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool2
+ kalman_fitter2.origin = origin
+ ckf.KalmanFitterTool2 = kalman_fitter2
+
+ ckf.TrackSeed = track_seed_tool
+ ckf.ActsLogging = "INFO"
+ ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool
+ ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool
+ ckf.PerformanceWriter = False
+ # ckf.PerformanceWriterTool = trajectory_performance_writer_tool
+ ckf.origin = origin
+
+ ckf.nMax = 10
+ ckf.chi2Max = 25
+ acc.addEventAlgo(ckf)
+ acc.merge(CKF2_OutputCfg(flags))
+ return acc
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 14d6b22aa987975dddc87875cc0dbb2fef7b7be2..f8f693b47c35aeb0cf938fa1932f6f384171b16b 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -277,9 +277,9 @@ StatusCode CKF2::execute() {
std::unique_ptr<Trk::Track> track = makeTrack(geoctx, traj);
if (track) {
// outputTracks->push_back(std::move(track));
- std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, *track, trajectories, Acts::BoundVector::Zero());
+ std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, *track, trajectories, Acts::BoundVector::Zero(), m_isMC);
if (track2) {
- std::unique_ptr<Trk::Track> track3 = m_kalmanFitterTool2->fit(ctx, *track2, trajectories, Acts::BoundVector::Zero());
+ std::unique_ptr<Trk::Track> track3 = m_kalmanFitterTool2->fit(ctx, *track2, trajectories, Acts::BoundVector::Zero(), m_isMC);
outputTracks->push_back(std::move(track3));
}
}
@@ -287,10 +287,10 @@ StatusCode CKF2::execute() {
// run the performance writer
if (m_statesWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, selectedTrajectories));
+ ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, selectedTrajectories, m_isMC));
}
if (m_summaryWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, selectedTrajectories));
+ ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, selectedTrajectories, m_isMC));
}
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(geoctx, selectedTrajectories));
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index a9cbb92919799f2cabf183bfcba4ffe3f611d917..5cc5aa20a80dbab1ac2120c48536b2652ee3ac1b 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -104,7 +104,10 @@ StatusCode CircleFitTrackSeedTool::run() {
// create seeds from four stations
// go(segments, combination, seeds, 0, 4);
// create seeds from three stations
- go(segments, combination, seeds, 0, 3);
+ go(segments, combination, seeds, 0, 4);
+ if (seeds.size() < 2) {
+ go(segments, combination, seeds, 0, 3);
+ }
// create seeds from two stations
if (seeds.size() < 2) {
go(segments, combination, seeds, 0, 2);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
index 2ee66e278cd3488a4c99611c11d73a21bac2b3d2..71afcc4431fabb6064af9601d31dd783de497a4f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
@@ -171,10 +171,10 @@ StatusCode CombinatorialKalmanFilterAlg::execute() {
// run the performance writer
if (m_statesWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, trajectories));
+ ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, trajectories, m_isMC));
}
if (m_summaryWriter && !m_noDiagnostics) {
- ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, trajectories));
+ ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, trajectories, m_isMC));
}
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(geoctx, trajectories));
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
index 4525f42cd575b0c3eee67893c09fca0dadeea3d6..6cfb07d143b8b49652ad54ee768a1af2b0671fae 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
@@ -12,7 +12,7 @@ GhostBusters::GhostBusters(const std::string &name, ISvcLocator *pSvcLocator)
StatusCode GhostBusters::initialize() {
ATH_CHECK(m_trackCollection.initialize());
- ATH_CHECK(m_simDataCollectionKey.initialize());
+ // ATH_CHECK(m_simDataCollectionKey.initialize());
ATH_CHECK(m_outputTrackCollection.initialize());
ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
@@ -23,7 +23,7 @@ StatusCode GhostBusters::initialize() {
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("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));
@@ -38,15 +38,15 @@ StatusCode GhostBusters::execute(const EventContext &ctx) const {
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<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);
+ auto segment = Segment(track, m_idHelper);
segments[segment.station()].push_back(segment);
}
@@ -77,7 +77,7 @@ StatusCode GhostBusters::execute(const EventContext &ctx) const {
m_chi2 = segment.chi2();
m_station = segment.station();
m_size = segment.size();
- m_majorityHits = segment.majorityHits();
+ // m_majorityHits = segment.majorityHits();
m_isGhost = segment.isGhost();
if (nGoodSegments >= 2 && segment.size() == 4) m_isGhost = true;
m_tree->Fill();
@@ -97,8 +97,7 @@ StatusCode GhostBusters::finalize() {
}
-GhostBusters::Segment::Segment(const Trk::Track *track, const TrackerSimDataCollection &simDataCollection,
- const FaserSCT_ID *idHelper) {
+GhostBusters::Segment::Segment(const Trk::Track *track, const FaserSCT_ID *idHelper) {
m_track = track;
m_position = track->trackParameters()->front()->position();
m_chi2 = track->fitQuality()->chiSquared();
@@ -111,5 +110,5 @@ GhostBusters::Segment::Segment(const Trk::Track *track, const TrackerSimDataColl
}
}
m_size = clusters.size();
- identifyContributingParticles(simDataCollection, clusters, m_particleHitCounts);
+ // identifyContributingParticles(simDataCollection, clusters, m_particleHitCounts);
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index a8e0ad41de7b668fd360265dcf5ebe5bd12178e6..4d696d19534865d8c5d62c32375b948e1352eb68 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -41,7 +41,7 @@ StatusCode KalmanFitterTool::finalize() {
std::unique_ptr<Trk::Track>
KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
std::vector<FaserActsRecMultiTrajectory> &trajectories,
- const Acts::BoundVector& inputVector = Acts::BoundVector::Zero()) const {
+ const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(), bool isMC=false) const {
std::unique_ptr<Trk::Track> newTrack = nullptr;
std::vector<FaserActsRecMultiTrajectory> myTrajectories;
@@ -108,10 +108,10 @@ KalmanFitterTool::fit(const EventContext &ctx, const Trk::Track &inputTrack,
}
if (m_statesWriter && !m_noDiagnostics) {
- StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories);
+ StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories, isMC);
}
if (m_summaryWriter && !m_noDiagnostics) {
- StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories);
+ StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories, isMC);
}
return newTrack;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
index 634e356a887d45afa6ddf7c06a1d8e4566d0ef8b..3964e0f0b98796ff840ee5597da6a8a75d0901c9 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
@@ -13,6 +13,9 @@
#include <TFile.h>
#include <TTree.h>
+constexpr float NaNfloat = std::numeric_limits<float>::quiet_NaN();
+constexpr float NaNint = std::numeric_limits<int>::quiet_NaN();
+
using Acts::VectorHelpers::eta;
using Acts::VectorHelpers::perp;
using Acts::VectorHelpers::phi;
@@ -194,6 +197,7 @@ StatusCode RootTrajectoryStatesWriterTool::initialize() {
m_outputTree->Branch("chi2", &m_chi2);
}
+
return StatusCode::SUCCESS;
}
@@ -207,51 +211,46 @@ StatusCode RootTrajectoryStatesWriterTool::finalize() {
return StatusCode::SUCCESS;
}
-StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories) const {
+StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const {
if (m_outputFile == nullptr)
return StatusCode::SUCCESS;
+ // Get the event number
const EventContext& ctx = Gaudi::Hive::currentContext();
+ m_eventNr = ctx.eventID().event_number();
- SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
- ATH_CHECK(mcEvents.isValid());
- if (mcEvents->size() != 1) {
- ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
- return StatusCode::FAILURE;
- }
- ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
+ std::vector<ParticleHitCount> particleHitCounts;
std::map<int, const HepMC::GenParticle*> particles {};
- for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
- particles[particle->barcode()] = particle;
- }
-
- SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
- ATH_CHECK(simData.isValid());
-
- SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
- ATH_CHECK(siHitCollection.isValid());
std::map<std::pair<int, Identifier>, const FaserSiHit*> siHitMap;
- for (const FaserSiHit& hit : *siHitCollection) {
- int barcode = hit.trackNumber();
- Identifier id = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
- siHitMap[std::make_pair(barcode, id)] = &hit;
- }
- // using HitParticlesMap = IndexMultimap<ActsFatras::Barcode>;
- // using HitSimHitsMap = IndexMultimap<Index>;
+ std::shared_ptr<TrackerSimDataCollection> simData {nullptr};
+
+ if (isMC) {
+ SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
+ ATH_CHECK(mcEvents.isValid());
+ if (mcEvents->size() != 1) {
+ ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
+ return StatusCode::FAILURE;
+ }
+ ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
+ for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ particles[particle->barcode()] = particle;
+ }
- // Read additional input collections
- // const auto& particles = ctx.eventStore.get<SimParticleContainer>(m_cfg.inputParticles);
- // const auto& simHits = ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimHits);
- // const auto& hitParticlesMap = ctx.eventStore.get<HitParticlesMap>(m_cfg.inputMeasurementParticlesMap);
- // const auto& hitSimHitsMap = ctx.eventStore.get<HitSimHitsMap>(m_cfg.inputMeasurementSimHitsMap);
+ SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
+ ATH_CHECK(simDataHandle.isValid());
+ simData = std::make_shared<TrackerSimDataCollection>(*simDataHandle);
- // For each particle within a track, how many hits did it contribute
- std::vector<ParticleHitCount> particleHitCounts;
+ SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
+ ATH_CHECK(siHitCollection.isValid());
+ for (const FaserSiHit& hit : *siHitCollection) {
+ int barcode = hit.trackNumber();
+ Identifier id = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
+ siHitMap[std::make_pair(barcode, id)] = &hit;
+ }
+ }
- // Get the event number
- m_eventNr = ctx.eventID().event_number();
// Loop over the trajectories
for (size_t itraj = 0; itraj < trajectories.size(); ++itraj) {
@@ -281,28 +280,38 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
m_nStates = trajState.nStates;
// Get the majority truth particle to this track
- int barcode;
- int truthQ = 1.;
- float truthMomentum = 1;
- identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
- if (not particleHitCounts.empty()) {
- // Get the barcode of the majority truth particle
- barcode = particleHitCounts.front().particleId;
- // Find the truth particle via the barcode
- auto ip = particles.find(barcode);
- if (ip != particles.end()) {
- const auto& particle = ip->second;
- ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
- // Get the truth particle charge
- // FIXME find better way to access charge of simulated particle, this does not work for
- // pions which have a positive pdg code (211) and positive charge
- truthQ = particle->pdg_id() > 0 ? -1 : 1;
- truthMomentum = particle->momentum().rho() * m_MeV2GeV;
- } else {
- ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
+ int barcode = NaNfloat;
+ int truthQ = NaNfloat;
+ float truthMomentum = NaNfloat;
+ float truthLOC0 = NaNfloat;
+ float truthLOC1 = NaNfloat;
+ float truthPHI = NaNfloat;
+ float truthTHETA = NaNfloat;
+ float truthQOP = NaNfloat;
+ float truthTIME = NaNfloat;
+
+ if (isMC) {
+ truthQ = 1;
+ truthMomentum = 1;
+ identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
+ if (not particleHitCounts.empty()) {
+ // Get the barcode of the majority truth particle
+ barcode = particleHitCounts.front().particleId;
+ // Find the truth particle via the barcode
+ auto ip = particles.find(barcode);
+ if (ip != particles.end()) {
+ const auto& particle = ip->second;
+ ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
+ // Get the truth particle charge
+ // FIXME find better way to access charge of simulated particle, this does not work for
+ // pions which have a positive pdg code (211) and positive charge
+ truthQ = particle->pdg_id() > 0 ? -1 : 1;
+ truthMomentum = particle->momentum().rho() * m_MeV2GeV;
+ } else {
+ ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
+ }
}
}
- float truthQOP = truthQ / truthMomentum;
// Get the trackStates on the trajectory
m_nParams = {0, 0, 0};
@@ -320,48 +329,63 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
// get the truth hits corresponding to this trackState
Identifier id = state.uncalibrated().hit()->identify();
Identifier waferId = m_idHelper->wafer_id(id);
- if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
- ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
- return true;
+ // if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
+ // ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
+ // return true;
+ // }
+
+ if (isMC) {
+ const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
+ HepGeom::Point3D localStartPos = siHit->localStartPosition();
+ HepGeom::Point3D localEndPos = siHit->localEndPosition();
+ HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
+ auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
+ const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
+ const HepGeom::Point3D<double> globalStartPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
+ const HepGeom::Point3D<double> globalEndPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
+ auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
+ auto globalDirection = globalEndPosition - globalStartPosition;
+ auto truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
+ auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
+ // FIXME get truthQOP for each state
+
+ // fill the truth hit info
+ m_t_x.push_back(truthPos[Acts::ePos0]);
+ m_t_y.push_back(truthPos[Acts::ePos1]);
+ m_t_z.push_back(truthPos[Acts::ePos2]);
+ m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
+ m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
+ m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
+
+ // get the truth track parameter at this track State
+ float truthLOC0 = truthLocal[Acts::ePos0];
+ float truthLOC1 = truthLocal[Acts::ePos1];
+ float truthTIME = siHit->meanTime();
+ float truthPHI = phi(truthUnitDir);
+ float truthTHETA = theta(truthUnitDir);
+
+ // fill the truth track parameter at this track State
+ m_t_eLOC0.push_back(truthLOC0);
+ m_t_eLOC1.push_back(truthLOC1);
+ m_t_ePHI.push_back(truthPHI);
+ m_t_eTHETA.push_back(truthTHETA);
+ m_t_eQOP.push_back(truthQOP);
+ m_t_eT.push_back(truthTIME);
+ } else {
+ m_t_x.push_back(NaNfloat);
+ m_t_y.push_back(NaNfloat);
+ m_t_z.push_back(NaNfloat);
+ m_t_dx.push_back(NaNfloat);
+ m_t_dy.push_back(NaNfloat);
+ m_t_dz.push_back(NaNfloat);
+ m_t_eLOC0.push_back(NaNfloat);
+ m_t_eLOC1.push_back(NaNfloat);
+ m_t_ePHI.push_back(NaNfloat);
+ m_t_eTHETA.push_back(NaNfloat);
+ m_t_eT.push_back(NaNfloat);
}
- const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
- HepGeom::Point3D localStartPos = siHit->localStartPosition();
- HepGeom::Point3D localEndPos = siHit->localEndPosition();
- HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
- auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
- const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
- const HepGeom::Point3D<double> globalStartPosition =
- Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
- const HepGeom::Point3D<double> globalEndPosition =
- Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
- auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
- auto globalDirection = globalEndPosition - globalStartPosition;
- auto truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
- auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
- // FIXME get truthQOP for each state
-
- // fill the truth hit info
- m_t_x.push_back(truthPos[Acts::ePos0]);
- m_t_y.push_back(truthPos[Acts::ePos1]);
- m_t_z.push_back(truthPos[Acts::ePos2]);
- m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
- m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
- m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
-
- // get the truth track parameter at this track State
- float truthLOC0 = truthLocal[Acts::ePos0];
- float truthLOC1 = truthLocal[Acts::ePos1];
- float truthTIME = siHit->meanTime();
- float truthPHI = phi(truthUnitDir);
- float truthTHETA = theta(truthUnitDir);
-
- // fill the truth track parameter at this track State
- m_t_eLOC0.push_back(truthLOC0);
- m_t_eLOC1.push_back(truthLOC1);
- m_t_ePHI.push_back(truthPHI);
- m_t_eTHETA.push_back(truthTHETA);
- m_t_eQOP.push_back(truthQOP);
- m_t_eT.push_back(truthTIME);
// get the geometry ID
auto geoID = surface.geometryId();
@@ -433,20 +457,20 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
auto res = state.effectiveCalibrated() - H * parameters;
m_res_x_hit.push_back(res[Acts::eBoundLoc0]);
// m_res_y_hit.push_back(res[Acts::eBoundLoc1]);
- m_res_y_hit.push_back(-99.);
+ m_res_y_hit.push_back(NaNfloat);
m_err_x_hit.push_back(
sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
// m_err_y_hit.push_back(
// sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_err_x_hit.push_back(-99.);
- m_res_y_hit.push_back(-99.);
+ m_err_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
m_pull_x_hit.push_back(
res[Acts::eBoundLoc0] /
sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
// m_pull_y_hit.push_back(
// res[Acts::eBoundLoc1] /
// sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_pull_y_hit.push_back(-99);
+ m_pull_y_hit.push_back(NaNfloat);
m_dim_hit.push_back(state.calibratedSize());
}
@@ -459,52 +483,94 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
m_eT[ipar].push_back(parameters[Acts::eBoundTime]);
// track parameters residual
- m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] -
- truthLOC0);
- m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] -
- truthLOC1);
- float resPhi = Acts::detail::difference_periodic<float>(
- parameters[Acts::eBoundPhi], truthPHI,
- static_cast<float>(2 * M_PI));
- m_res_ePHI[ipar].push_back(resPhi);
- m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] -
- truthTHETA);
- m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] -
- truthQOP);
- m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
-
- // track parameters error
- m_err_eLOC0[ipar].push_back(
- sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
- m_err_eLOC1[ipar].push_back(
- sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_err_ePHI[ipar].push_back(
- sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
- m_err_eTHETA[ipar].push_back(
- sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
- m_err_eQOP[ipar].push_back(
- sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
- m_err_eT[ipar].push_back(
- sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
-
- // track parameters pull
- m_pull_eLOC0[ipar].push_back(
- (parameters[Acts::eBoundLoc0] - truthLOC0) /
- sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
- m_pull_eLOC1[ipar].push_back(
- (parameters[Acts::eBoundLoc1] - truthLOC1) /
- sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
- m_pull_ePHI[ipar].push_back(
- resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
- m_pull_eTHETA[ipar].push_back(
- (parameters[Acts::eBoundTheta] - truthTHETA) /
- sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
- m_pull_eQOP[ipar].push_back(
- (parameters[Acts::eBoundQOverP] - truthQOP) /
- sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
- m_pull_eT[ipar].push_back(
- (parameters[Acts::eBoundTime] - truthTIME) /
- sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+ float resPhi;
+ if (isMC) {
+ m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
+ m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
+ resPhi = Acts::detail::difference_periodic<float>( parameters[Acts::eBoundPhi], truthPHI,
+ static_cast<float>(2 * M_PI));
+ m_res_ePHI[ipar].push_back(resPhi);
+ m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] - truthTHETA);
+ m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] - truthQOP);
+ m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
+
+ // track parameters error
+ m_err_eLOC0[ipar].push_back(
+ sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ m_err_eLOC1[ipar].push_back(
+ sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+ m_err_ePHI[ipar].push_back(
+ sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+ m_err_eTHETA[ipar].push_back(
+ sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+ m_err_eQOP[ipar].push_back(
+ sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+ m_err_eT[ipar].push_back(
+ sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+
+ // track parameters pull
+ m_pull_eLOC0[ipar].push_back(
+ (parameters[Acts::eBoundLoc0] - truthLOC0) /
+ sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ m_pull_eLOC1[ipar].push_back(
+ (parameters[Acts::eBoundLoc1] - truthLOC1) /
+ sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+ m_pull_ePHI[ipar].push_back(
+ resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+ m_pull_eTHETA[ipar].push_back(
+ (parameters[Acts::eBoundTheta] - truthTHETA) /
+ sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+ m_pull_eQOP[ipar].push_back(
+ (parameters[Acts::eBoundQOverP] - truthQOP) /
+ sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+ m_pull_eT[ipar].push_back(
+ (parameters[Acts::eBoundTime] - truthTIME) /
+ sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+ } else {
+ if (ipar == 0) {
+ // push default values if no track parameters
+ m_res_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
+ m_err_x_hit.push_back(NaNfloat);
+ m_err_y_hit.push_back(NaNfloat);
+ m_pull_x_hit.push_back(NaNfloat);
+ m_pull_y_hit.push_back(NaNfloat);
+ m_dim_hit.push_back(NaNfloat);
+ }
+ // push default values if no track parameters
+ // m_eLOC0[ipar].push_back(NaNfloat);
+ // m_eLOC1[ipar].push_back(NaNfloat);
+ // m_ePHI[ipar].push_back(NaNfloat);
+ // m_eTHETA[ipar].push_back(NaNfloat);
+ // m_eQOP[ipar].push_back(NaNfloat);
+ // m_eT[ipar].push_back(NaNfloat);
+ m_res_eLOC0[ipar].push_back(NaNfloat);
+ m_res_eLOC1[ipar].push_back(NaNfloat);
+ m_res_ePHI[ipar].push_back(NaNfloat);
+ m_res_eTHETA[ipar].push_back(NaNfloat);
+ m_res_eQOP[ipar].push_back(NaNfloat);
+ m_res_eT[ipar].push_back(NaNfloat);
+ m_err_eLOC0[ipar].push_back(NaNfloat);
+ m_err_eLOC1[ipar].push_back(NaNfloat);
+ m_err_ePHI[ipar].push_back(NaNfloat);
+ m_err_eTHETA[ipar].push_back(NaNfloat);
+ m_err_eQOP[ipar].push_back(NaNfloat);
+ m_err_eT[ipar].push_back(NaNfloat);
+ m_pull_eLOC0[ipar].push_back(NaNfloat);
+ m_pull_eLOC1[ipar].push_back(NaNfloat);
+ m_pull_ePHI[ipar].push_back(NaNfloat);
+ m_pull_eTHETA[ipar].push_back(NaNfloat);
+ m_pull_eQOP[ipar].push_back(NaNfloat);
+ m_pull_eT[ipar].push_back(NaNfloat);
+ // m_x[ipar].push_back(NaNfloat);
+ // m_y[ipar].push_back(NaNfloat);
+ // m_z[ipar].push_back(NaNfloat);
+ // m_px[ipar].push_back(NaNfloat);
+ // m_py[ipar].push_back(NaNfloat);
+ // m_pz[ipar].push_back(NaNfloat);
+ // m_pT[ipar].push_back(NaNfloat);
+ // m_eta[ipar].push_back(NaNfloat);
+ }
// further track parameter info
Acts::FreeVector freeParams =
@@ -524,47 +590,47 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
} else {
if (ipar == 0) {
// push default values if no track parameters
- m_res_x_hit.push_back(-99.);
- m_res_y_hit.push_back(-99.);
- m_err_x_hit.push_back(-99.);
- m_err_y_hit.push_back(-99.);
- m_pull_x_hit.push_back(-99.);
- m_pull_y_hit.push_back(-99.);
- m_dim_hit.push_back(-99.);
+ m_res_x_hit.push_back(NaNfloat);
+ m_res_y_hit.push_back(NaNfloat);
+ m_err_x_hit.push_back(NaNfloat);
+ m_err_y_hit.push_back(NaNfloat);
+ m_pull_x_hit.push_back(NaNfloat);
+ m_pull_y_hit.push_back(NaNfloat);
+ m_dim_hit.push_back(NaNfloat);
}
// push default values if no track parameters
- m_eLOC0[ipar].push_back(-99.);
- m_eLOC1[ipar].push_back(-99.);
- m_ePHI[ipar].push_back(-99.);
- m_eTHETA[ipar].push_back(-99.);
- m_eQOP[ipar].push_back(-99.);
- m_eT[ipar].push_back(-99.);
- m_res_eLOC0[ipar].push_back(-99.);
- m_res_eLOC1[ipar].push_back(-99.);
- m_res_ePHI[ipar].push_back(-99.);
- m_res_eTHETA[ipar].push_back(-99.);
- m_res_eQOP[ipar].push_back(-99.);
- m_res_eT[ipar].push_back(-99.);
- m_err_eLOC0[ipar].push_back(-99);
- m_err_eLOC1[ipar].push_back(-99);
- m_err_ePHI[ipar].push_back(-99);
- m_err_eTHETA[ipar].push_back(-99);
- m_err_eQOP[ipar].push_back(-99);
- m_err_eT[ipar].push_back(-99);
- m_pull_eLOC0[ipar].push_back(-99.);
- m_pull_eLOC1[ipar].push_back(-99.);
- m_pull_ePHI[ipar].push_back(-99.);
- m_pull_eTHETA[ipar].push_back(-99.);
- m_pull_eQOP[ipar].push_back(-99.);
- m_pull_eT[ipar].push_back(-99.);
- m_x[ipar].push_back(-99.);
- m_y[ipar].push_back(-99.);
- m_z[ipar].push_back(-99.);
- m_px[ipar].push_back(-99.);
- m_py[ipar].push_back(-99.);
- m_pz[ipar].push_back(-99.);
- m_pT[ipar].push_back(-99.);
- m_eta[ipar].push_back(-99.);
+ m_eLOC0[ipar].push_back(NaNfloat);
+ m_eLOC1[ipar].push_back(NaNfloat);
+ m_ePHI[ipar].push_back(NaNfloat);
+ m_eTHETA[ipar].push_back(NaNfloat);
+ m_eQOP[ipar].push_back(NaNfloat);
+ m_eT[ipar].push_back(NaNfloat);
+ m_res_eLOC0[ipar].push_back(NaNfloat);
+ m_res_eLOC1[ipar].push_back(NaNfloat);
+ m_res_ePHI[ipar].push_back(NaNfloat);
+ m_res_eTHETA[ipar].push_back(NaNfloat);
+ m_res_eQOP[ipar].push_back(NaNfloat);
+ m_res_eT[ipar].push_back(NaNfloat);
+ m_err_eLOC0[ipar].push_back(NaNfloat);
+ m_err_eLOC1[ipar].push_back(NaNfloat);
+ m_err_ePHI[ipar].push_back(NaNfloat);
+ m_err_eTHETA[ipar].push_back(NaNfloat);
+ m_err_eQOP[ipar].push_back(NaNfloat);
+ m_err_eT[ipar].push_back(NaNfloat);
+ m_pull_eLOC0[ipar].push_back(NaNfloat);
+ m_pull_eLOC1[ipar].push_back(NaNfloat);
+ m_pull_ePHI[ipar].push_back(NaNfloat);
+ m_pull_eTHETA[ipar].push_back(NaNfloat);
+ m_pull_eQOP[ipar].push_back(NaNfloat);
+ m_pull_eT[ipar].push_back(NaNfloat);
+ m_x[ipar].push_back(NaNfloat);
+ m_y[ipar].push_back(NaNfloat);
+ m_z[ipar].push_back(NaNfloat);
+ m_px[ipar].push_back(NaNfloat);
+ m_py[ipar].push_back(NaNfloat);
+ m_pz[ipar].push_back(NaNfloat);
+ m_pT[ipar].push_back(NaNfloat);
+ m_eta[ipar].push_back(NaNfloat);
}
// fill the track parameters status
m_hasParams[ipar].push_back(hasParams[ipar]);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
index bb916b9bfefeaf47c013b29122c12eb48f17abf7..b4b97c1586048f027a7d7e353e4da88b8df2cfb9 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
@@ -15,7 +15,6 @@
#include <TTree.h>
/// NaN values for TTree variables
-constexpr double NaNdouble = std::numeric_limits<double>::quiet_NaN();
constexpr float NaNfloat = std::numeric_limits<float>::quiet_NaN();
constexpr float NaNint = std::numeric_limits<int>::quiet_NaN();
@@ -119,32 +118,35 @@ StatusCode RootTrajectorySummaryWriterTool::initialize() {
StatusCode RootTrajectorySummaryWriterTool::finalize() {
- if (!m_noDiagnostics) {
- m_outputFile->cd();
- m_outputTree->Write();
- m_outputFile->Close();
- }
+ m_outputFile->cd();
+ m_outputTree->Write();
+ m_outputFile->Close();
return StatusCode::SUCCESS;
}
StatusCode RootTrajectorySummaryWriterTool::write(
- const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories) const {
+ const Acts::GeometryContext& geoContext, const TrajectoriesContainer& trajectories, bool isMC) const {
EventContext ctx = Gaudi::Hive::currentContext();
- SG::ReadHandle<TrackerSimDataCollection> simData {m_simDataCollectionKey, ctx};
- ATH_CHECK(simData.isValid());
+ std::shared_ptr<TrackerSimDataCollection> simData {nullptr};
+ std::map<int, const HepMC::GenParticle*> particles {};
- SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
- ATH_CHECK(mcEvents.isValid());
- if (mcEvents->size() != 1) {
- ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
- return StatusCode::FAILURE;
- }
+ if (isMC) {
+ SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
+ ATH_CHECK(simDataHandle.isValid());
+ simData = std::make_shared<TrackerSimDataCollection>(*simDataHandle);
- std::map<int, const HepMC::GenParticle*> particles {};
- for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
- particles[particle->barcode()] = particle;
+ SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
+ ATH_CHECK(mcEvents.isValid());
+ if (mcEvents->size() != 1) {
+ ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
+ return StatusCode::FAILURE;
+ }
+
+ for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
+ particles[particle->barcode()] = particle;
+ }
}
// For each particle within a track, how many hits did it contribute
@@ -219,69 +221,71 @@ StatusCode RootTrajectorySummaryWriterTool::write(
pSurface = const_cast<Acts::Surface*>(&boundParam.referenceSurface());
}
- // Get the majority truth particle to this track
- ATH_MSG_VERBOSE("get majority truth particle");
- identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
- for (const auto& particle : particleHitCounts) {
- ATH_MSG_VERBOSE(particle.particleId << ": " << particle.hitCount << " hits");
- }
+ if (isMC) {
+ // Get the majority truth particle to this track
+ ATH_MSG_VERBOSE("get majority truth particle");
+ identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
+ for (const auto& particle : particleHitCounts) {
+ ATH_MSG_VERBOSE(particle.particleId << ": " << particle.hitCount << " hits");
+ }
- bool foundMajorityParticle = false;
- // Get the truth particle info
- if (not particleHitCounts.empty()) {
- // Get the barcode of the majority truth particle
- majorityParticleId = particleHitCounts.front().particleId;
- nMajorityHits = particleHitCounts.front().hitCount;
-
- // Find the truth particle via the barcode
- auto ip = particles.find(majorityParticleId);
- if (ip != particles.end()) {
- foundMajorityParticle = true;
-
- const HepMC::GenParticle* particle = ip->second;
- ATH_MSG_DEBUG("Find the truth particle with barcode = " << majorityParticleId);
-
- // extrapolate parameters from vertex to reference surface at origin.
- std::unique_ptr<const Acts::BoundTrackParameters> truthParameters
- = extrapolateToReferenceSurface(ctx, particle);
- if (!truthParameters) {
- continue;
- }
- // Get the truth particle info at vertex
- const HepMC::GenVertex* vertex = particle->production_vertex();
- t_p = truthParameters->momentum().mag();
- t_charge = truthParameters->charge();
- t_time = truthParameters->time();
- t_vx = truthParameters->position(geoContext).x();
- t_vy = truthParameters->position(geoContext).y();
- t_vz = truthParameters->position(geoContext).z();
- t_px = truthParameters->momentum().x();
- t_py = truthParameters->momentum().y();
- t_pz = truthParameters->momentum().z();
- t_theta = theta(truthParameters->momentum().normalized());
- t_phi = phi(truthParameters->momentum().normalized());
- t_eta = eta(truthParameters->momentum().normalized());
- t_pT = t_p * perp(truthParameters->momentum().normalized());
-
- auto unitDirection = Acts::Vector3(t_px, t_py, t_pz).normalized();
- auto positon = Acts::Vector3 (t_vx, t_vy, t_vz);
-
- if (pSurface) {
- // get the truth perigee parameter
- auto lpResult = pSurface->globalToLocal(geoContext, positon, unitDirection);
- if (lpResult.ok()) {
- t_vlx = lpResult.value()[Acts::BoundIndices::eBoundLoc0];
- t_vly = lpResult.value()[Acts::BoundIndices::eBoundLoc1];
- } else {
- ATH_MSG_ERROR("Global to local transformation did not succeed.");
+ bool foundMajorityParticle = false;
+ // Get the truth particle info
+ if (not particleHitCounts.empty()) {
+ // Get the barcode of the majority truth particle
+ majorityParticleId = particleHitCounts.front().particleId;
+ nMajorityHits = particleHitCounts.front().hitCount;
+
+ // Find the truth particle via the barcode
+ auto ip = particles.find(majorityParticleId);
+ if (ip != particles.end()) {
+ foundMajorityParticle = true;
+
+ const HepMC::GenParticle* particle = ip->second;
+ ATH_MSG_DEBUG("Find the truth particle with barcode = " << majorityParticleId);
+
+ // extrapolate parameters from vertex to reference surface at origin.
+ std::unique_ptr<const Acts::BoundTrackParameters> truthParameters
+ = extrapolateToReferenceSurface(ctx, particle);
+ if (!truthParameters) {
+ continue;
+ }
+ // Get the truth particle info at vertex
+ const HepMC::GenVertex* vertex = particle->production_vertex();
+ t_p = truthParameters->momentum().mag();
+ t_charge = truthParameters->charge();
+ t_time = truthParameters->time();
+ t_vx = truthParameters->position(geoContext).x();
+ t_vy = truthParameters->position(geoContext).y();
+ t_vz = truthParameters->position(geoContext).z();
+ t_px = truthParameters->momentum().x();
+ t_py = truthParameters->momentum().y();
+ t_pz = truthParameters->momentum().z();
+ t_theta = theta(truthParameters->momentum().normalized());
+ t_phi = phi(truthParameters->momentum().normalized());
+ t_eta = eta(truthParameters->momentum().normalized());
+ t_pT = t_p * perp(truthParameters->momentum().normalized());
+
+ auto unitDirection = Acts::Vector3(t_px, t_py, t_pz).normalized();
+ auto positon = Acts::Vector3 (t_vx, t_vy, t_vz);
+
+ if (pSurface) {
+ // get the truth perigee parameter
+ auto lpResult = pSurface->globalToLocal(geoContext, positon, unitDirection);
+ if (lpResult.ok()) {
+ t_vlx = lpResult.value()[Acts::BoundIndices::eBoundLoc0];
+ t_vly = lpResult.value()[Acts::BoundIndices::eBoundLoc1];
+ } else {
+ ATH_MSG_ERROR("Global to local transformation did not succeed.");
+ }
}
+ } else {
+ ATH_MSG_WARNING("Truth particle with barcode = " << majorityParticleId << " not found in the input collection!");
}
- } else {
- ATH_MSG_WARNING("Truth particle with barcode = " << majorityParticleId << " not found in the input collection!");
}
- }
- if (not foundMajorityParticle) {
- ATH_MSG_WARNING("Truth particle for mj " << itraj << " subtraj " << isubtraj << " not found!");
+ if (not foundMajorityParticle) {
+ ATH_MSG_WARNING("Truth particle for mj " << itraj << " subtraj " << isubtraj << " not found!");
+ }
}
// Push the corresponding truth particle info for the track.
diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py b/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py
new file mode 100644
index 0000000000000000000000000000000000000000..a22a4c407ffde297502d71b1d0a266489e6832c0
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/test/TI12CKF2.py
@@ -0,0 +1,62 @@
+#!/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 AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
+from FaserActsKalmanFilter.TI12CKF2Config import TI12CKF2Cfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['threeStationRun6833.raw']
+ConfigFlags.Output.ESDFileName = "CKF.ESD.pool.root"
+ConfigFlags.addFlag("Output.xAODFileName", f"CKF.xAOD.root")
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
+ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"
+ConfigFlags.Input.ProjectName = "data21"
+ConfigFlags.Input.isMC = False
+ConfigFlags.GeoModel.FaserVersion = "FASER-02"
+ConfigFlags.Common.isOnline = False
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.Detector.GeometryFaserSCT = True
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolWriteCfg(ConfigFlags))
+acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_LEVELMODE_RDOs", ClusterToolTimingPattern="X1X"))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+acc.merge(GhostBustersCfg(ConfigFlags))
+acc.merge(TI12CKF2Cfg(ConfigFlags, noDiagnostics=False))
+acc.getEventAlgo("CKF2").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()
+
+# Hack to avoid problem with our use of MC databases when isMC = False
+replicaSvc = acc.getService("DBReplicaSvc")
+replicaSvc.COOLSQLiteVetoPattern = ""
+replicaSvc.UseCOOLSQLite = True
+replicaSvc.UseCOOLFrontier = False
+replicaSvc.UseGeomSQLite = True
+
+sc = acc.run(maxEvents=20)
+sys.exit(not sc.isSuccess())