diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d3466a6e166471480b9d7159868e4693adc79bb
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h
@@ -0,0 +1,111 @@
+#ifndef FASERACTSKALMANFILTER_ROOTTRAJECTORYSTATESWRITERTOOL_H
+#define FASERACTSKALMANFILTER_ROOTTRAJECTORYSTATESWRITERTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include <array>
+#include <string>
+#include <vector>
+
+
+class TFile;
+class TTree;
+struct FaserActsRecMultiTrajectory;
+using TrajectoriesContainer = std::vector<FaserActsRecMultiTrajectory>;
+
+class RootTrajectoryStatesWriterTool : public AthAlgTool {
+public:
+ RootTrajectoryStatesWriterTool(const std::string& type, const std::string& name,
+ const IInterface* parent);
+ ~RootTrajectoryStatesWriterTool() override = default;
+
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+ StatusCode write(TrajectoriesContainer trajectories, Acts::GeometryContext geoContext,
+ std::vector<Acts::CurvilinearTrackParameters> traj) const;
+
+private:
+ Gaudi::Property<std::string> m_filePath{this, "FilePath", "Trajectories.root", "Output root file"};
+ Gaudi::Property<std::string> m_treeName{this, "TreeName", "tree", "Tree name"};
+ TFile* m_file;
+ TTree* m_tree;
+
+ mutable uint32_t m_eventNr{0}; ///< the event number
+ mutable uint32_t m_multiTrajNr{0}; ///< the multi-trajectory number
+ mutable unsigned int m_subTrajNr{0}; ///< the multi-trajectory sub-trajectory number
+
+ mutable std::vector<float> m_t_x; ///< Global truth hit position x
+ mutable std::vector<float> m_t_y; ///< Global truth hit position y
+ mutable std::vector<float> m_t_z; ///< Global truth hit position z
+ mutable std::vector<float> m_t_r; ///< Global truth hit position r
+ mutable std::vector<float> m_t_dx; ///< Truth particle direction x at global hit position
+ mutable std::vector<float> m_t_dy; ///< Truth particle direction y at global hit position
+ mutable std::vector<float> m_t_dz; ///< Truth particle direction z at global hit position
+
+ mutable std::vector<float> m_t_eLOC0; ///< truth parameter eBoundLoc0
+ mutable std::vector<float> m_t_eLOC1; ///< truth parameter eBoundLoc1
+ mutable std::vector<float> m_t_ePHI; ///< truth parameter ePHI
+ mutable std::vector<float> m_t_eTHETA; ///< truth parameter eTHETA
+ mutable std::vector<float> m_t_eQOP; ///< truth parameter eQOP
+ mutable std::vector<float> m_t_eT; ///< truth parameter eT
+
+ mutable unsigned int m_nStates{0}; ///< number of all states
+ mutable unsigned int m_nMeasurements{0}; ///< number of states with measurements
+ mutable std::vector<int> m_volumeID; ///< volume identifier
+ mutable std::vector<int> m_layerID; ///< layer identifier
+ mutable std::vector<int> m_moduleID; ///< surface identifier
+ mutable std::vector<float> m_pathLength; ///< path length
+ mutable std::vector<float> m_lx_hit; ///< uncalibrated measurement local x
+ mutable std::vector<float> m_ly_hit; ///< uncalibrated measurement local y
+ mutable std::vector<float> m_x_hit; ///< uncalibrated measurement global x
+ mutable std::vector<float> m_y_hit; ///< uncalibrated measurement global y
+ mutable std::vector<float> m_z_hit; ///< uncalibrated measurement global z
+ mutable std::vector<float> m_res_x_hit; ///< hit residual x
+ mutable std::vector<float> m_res_y_hit; ///< hit residual y
+ mutable std::vector<float> m_err_x_hit; ///< hit err x
+ mutable std::vector<float> m_err_y_hit; ///< hit err y
+ mutable std::vector<float> m_pull_x_hit; ///< hit pull x
+ mutable std::vector<float> m_pull_y_hit; ///< hit pull y
+ mutable std::vector<int> m_dim_hit; ///< dimension of measurement
+
+ mutable std::array<int, 3> m_nParams; ///< number of states which have filtered/predicted/smoothed parameters
+ mutable std::array<std::vector<bool>, 3> m_hasParams; ///< status of the filtered/predicted/smoothed parameters
+ mutable std::array<std::vector<float>, 3> m_eLOC0; ///< predicted/filtered/smoothed parameter eLOC0
+ mutable std::array<std::vector<float>, 3> m_eLOC1; ///< predicted/filtered/smoothed parameter eLOC1
+ mutable std::array<std::vector<float>, 3> m_ePHI; ///< predicted/filtered/smoothed parameter ePHI
+ mutable std::array<std::vector<float>, 3> m_eTHETA; ///< predicted/filtered/smoothed parameter eTHETA
+ mutable std::array<std::vector<float>, 3> m_eQOP; ///< predicted/filtered/smoothed parameter eQOP
+ mutable std::array<std::vector<float>, 3> m_eT; ///< predicted/filtered/smoothed parameter eT
+ mutable std::array<std::vector<float>, 3> m_res_eLOC0; ///< predicted/filtered/smoothed parameter eLOC0 residual
+ mutable std::array<std::vector<float>, 3> m_res_eLOC1; ///< predicted/filtered/smoothed parameter eLOC1 residual
+ mutable std::array<std::vector<float>, 3> m_res_ePHI; ///< predicted/filtered/smoothed parameter ePHI residual
+ mutable std::array<std::vector<float>, 3> m_res_eTHETA; ///< predicted/filtered/smoothed parameter eTHETA residual
+ mutable std::array<std::vector<float>, 3> m_res_eQOP; ///< predicted/filtered/smoothed parameter eQOP residual
+ mutable std::array<std::vector<float>, 3> m_res_eT; ///< predicted/filtered/smoothed parameter eT residual
+ mutable std::array<std::vector<float>, 3> m_err_eLOC0; ///< predicted/filtered/smoothed parameter eLOC0 error
+ mutable std::array<std::vector<float>, 3> m_err_eLOC1; ///< predicted/filtered/smoothed parameter eLOC1 error
+ mutable std::array<std::vector<float>, 3> m_err_ePHI; ///< predicted/filtered/smoothed parameter ePHI error
+ mutable std::array<std::vector<float>, 3> m_err_eTHETA; ///< predicted/filtered/smoothed parameter eTHETA error
+ mutable std::array<std::vector<float>, 3> m_err_eQOP; ///< predicted/filtered/smoothed parameter eQOP error
+ mutable std::array<std::vector<float>, 3> m_err_eT; ///< predicted/filtered/smoothed parameter eT error
+ mutable std::array<std::vector<float>, 3> m_pull_eLOC0; ///< predicted/filtered/smoothed parameter eLOC0 pull
+ mutable std::array<std::vector<float>, 3> m_pull_eLOC1; ///< predicted/filtered/smoothed parameter eLOC1 pull
+ mutable std::array<std::vector<float>, 3> m_pull_ePHI; ///< predicted/filtered/smoothed parameter ePHI pull
+ mutable std::array<std::vector<float>, 3> m_pull_eTHETA; ///< predicted/filtered/smoothed parameter eTHETA pull
+ mutable std::array<std::vector<float>, 3> m_pull_eQOP; ///< predicted/filtered/smoothed parameter eQOP pull
+ mutable std::array<std::vector<float>, 3> m_pull_eT; ///< predicted/filtered/smoothed parameter eT pull
+ mutable std::array<std::vector<float>, 3> m_x; ///< predicted/filtered/smoothed parameter global x
+ mutable std::array<std::vector<float>, 3> m_y; ///< predicted/filtered/smoothed parameter global y
+ mutable std::array<std::vector<float>, 3> m_z; ///< predicted/filtered/smoothed parameter global z
+ mutable std::array<std::vector<float>, 3> m_px; ///< predicted/filtered/smoothed parameter px
+ mutable std::array<std::vector<float>, 3> m_py; ///< predicted/filtered/smoothed parameter py
+ mutable std::array<std::vector<float>, 3> m_pz; ///< predicted/filtered/smoothed parameter pz
+ mutable std::array<std::vector<float>, 3> m_eta; ///< predicted/filtered/smoothed parameter eta
+ mutable std::array<std::vector<float>, 3> m_pT; ///< predicted/filtered/smoothed parameter pT
+
+ mutable std::vector<float> m_chi2; ///< chisq from filtering
+};
+
+#endif // FASERACTSKALMANFILTER_ROOTTRAJECTORYSTATESWRITERTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ec5a075c651f1bead4dd623ba63bf710dffccc5e
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
@@ -0,0 +1,576 @@
+#include "FaserActsKalmanFilter/RootTrajectoryStatesWriterTool.h"
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/MultiTrajectoryHelpers.hpp"
+#include "Acts/EventData/detail/TransformationBoundToFree.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "FaserActsKalmanFilter/FaserActsRecMultiTrajectory.h"
+#include <TFile.h>
+#include <TTree.h>
+
+using Acts::VectorHelpers::eta;
+using Acts::VectorHelpers::perp;
+using Acts::VectorHelpers::phi;
+using Acts::VectorHelpers::theta;
+
+RootTrajectoryStatesWriterTool::RootTrajectoryStatesWriterTool(
+ const std::string& type, const std::string& name, const IInterface* parent)
+ : AthAlgTool(type, name, parent) {}
+
+StatusCode RootTrajectoryStatesWriterTool::initialize() {
+ std::string filePath = m_filePath;
+ std::string treeName = m_treeName;
+ m_file = TFile::Open(filePath.c_str(), "RECREATE");
+ if (m_file == nullptr) {
+ ATH_MSG_ERROR("Unable to open output file at " << m_filePath);
+ return StatusCode::FAILURE;
+ }
+ m_file->cd();
+ m_tree = new TTree(treeName.c_str(), treeName.c_str());
+ if (m_tree == nullptr) {
+ ATH_MSG_ERROR("Unable to create TTree");
+ return StatusCode::FAILURE;
+ }
+
+ m_tree = new TTree("tree", "tree");
+
+ m_tree->Branch("event_nr", &m_eventNr);
+ m_tree->Branch("multiTraj_nr", &m_multiTrajNr);
+ m_tree->Branch("subTraj_nr", &m_subTrajNr);
+
+ m_tree->Branch("t_x", &m_t_x);
+ m_tree->Branch("t_y", &m_t_y);
+ m_tree->Branch("t_z", &m_t_z);
+ m_tree->Branch("t_r", &m_t_r);
+ m_tree->Branch("t_dx", &m_t_dx);
+ m_tree->Branch("t_dy", &m_t_dy);
+ m_tree->Branch("t_dz", &m_t_dz);
+ m_tree->Branch("t_eLOC0", &m_t_eLOC0);
+ m_tree->Branch("t_eLOC1", &m_t_eLOC1);
+ m_tree->Branch("t_ePHI", &m_t_ePHI);
+ m_tree->Branch("t_eTHETA", &m_t_eTHETA);
+ m_tree->Branch("t_eQOP", &m_t_eQOP);
+ m_tree->Branch("t_eT", &m_t_eT);
+
+ m_tree->Branch("nStates", &m_nStates);
+ m_tree->Branch("nMeasurements", &m_nMeasurements);
+ m_tree->Branch("volume_id", &m_volumeID);
+ m_tree->Branch("layer_id", &m_layerID);
+ m_tree->Branch("module_id", &m_moduleID);
+ m_tree->Branch("pathLength", &m_pathLength);
+ m_tree->Branch("l_x_hit", &m_lx_hit);
+ m_tree->Branch("l_y_hit", &m_ly_hit);
+ m_tree->Branch("g_x_hit", &m_x_hit);
+ m_tree->Branch("g_y_hit", &m_y_hit);
+ m_tree->Branch("g_z_hit", &m_z_hit);
+ m_tree->Branch("res_x_hit", &m_res_x_hit);
+ m_tree->Branch("res_y_hit", &m_res_y_hit);
+ m_tree->Branch("err_x_hit", &m_err_x_hit);
+ m_tree->Branch("err_y_hit", &m_err_y_hit);
+ m_tree->Branch("pull_x_hit", &m_pull_x_hit);
+ m_tree->Branch("pull_y_hit", &m_pull_y_hit);
+ m_tree->Branch("dim_hit", &m_dim_hit);
+
+ m_tree->Branch("nPredicted", &m_nParams[0]);
+ m_tree->Branch("predicted", &m_hasParams[0]);
+ m_tree->Branch("eLOC0_prt", &m_eLOC0[0]);
+ m_tree->Branch("eLOC1_prt", &m_eLOC1[0]);
+ m_tree->Branch("ePHI_prt", &m_ePHI[0]);
+ m_tree->Branch("eTHETA_prt", &m_eTHETA[0]);
+ m_tree->Branch("eQOP_prt", &m_eQOP[0]);
+ m_tree->Branch("eT_prt", &m_eT[0]);
+ m_tree->Branch("res_eLOC0_prt", &m_res_eLOC0[0]);
+ m_tree->Branch("res_eLOC1_prt", &m_res_eLOC1[0]);
+ m_tree->Branch("res_ePHI_prt", &m_res_ePHI[0]);
+ m_tree->Branch("res_eTHETA_prt", &m_res_eTHETA[0]);
+ m_tree->Branch("res_eQOP_prt", &m_res_eQOP[0]);
+ m_tree->Branch("res_eT_prt", &m_res_eT[0]);
+ m_tree->Branch("err_eLOC0_prt", &m_err_eLOC0[0]);
+ m_tree->Branch("err_eLOC1_prt", &m_err_eLOC1[0]);
+ m_tree->Branch("err_ePHI_prt", &m_err_ePHI[0]);
+ m_tree->Branch("err_eTHETA_prt", &m_err_eTHETA[0]);
+ m_tree->Branch("err_eQOP_prt", &m_err_eQOP[0]);
+ m_tree->Branch("err_eT_prt", &m_err_eT[0]);
+ m_tree->Branch("pull_eLOC0_prt", &m_pull_eLOC0[0]);
+ m_tree->Branch("pull_eLOC1_prt", &m_pull_eLOC1[0]);
+ m_tree->Branch("pull_ePHI_prt", &m_pull_ePHI[0]);
+ m_tree->Branch("pull_eTHETA_prt", &m_pull_eTHETA[0]);
+ m_tree->Branch("pull_eQOP_prt", &m_pull_eQOP[0]);
+ m_tree->Branch("pull_eT_prt", &m_pull_eT[0]);
+ m_tree->Branch("g_x_prt", &m_x[0]);
+ m_tree->Branch("g_y_prt", &m_y[0]);
+ m_tree->Branch("g_z_prt", &m_z[0]);
+ m_tree->Branch("px_prt", &m_px[0]);
+ m_tree->Branch("py_prt", &m_py[0]);
+ m_tree->Branch("pz_prt", &m_pz[0]);
+ m_tree->Branch("eta_prt", &m_eta[0]);
+ m_tree->Branch("pT_prt", &m_pT[0]);
+
+ m_tree->Branch("nFiltered", &m_nParams[1]);
+ m_tree->Branch("filtered", &m_hasParams[1]);
+ m_tree->Branch("eLOC0_flt", &m_eLOC0[1]);
+ m_tree->Branch("eLOC1_flt", &m_eLOC1[1]);
+ m_tree->Branch("ePHI_flt", &m_ePHI[1]);
+ m_tree->Branch("eTHETA_flt", &m_eTHETA[1]);
+ m_tree->Branch("eQOP_flt", &m_eQOP[1]);
+ m_tree->Branch("eT_flt", &m_eT[1]);
+ m_tree->Branch("res_eLOC0_flt", &m_res_eLOC0[1]);
+ m_tree->Branch("res_eLOC1_flt", &m_res_eLOC1[1]);
+ m_tree->Branch("res_ePHI_flt", &m_res_ePHI[1]);
+ m_tree->Branch("res_eTHETA_flt", &m_res_eTHETA[1]);
+ m_tree->Branch("res_eQOP_flt", &m_res_eQOP[1]);
+ m_tree->Branch("res_eT_flt", &m_res_eT[1]);
+ m_tree->Branch("err_eLOC0_flt", &m_err_eLOC0[1]);
+ m_tree->Branch("err_eLOC1_flt", &m_err_eLOC1[1]);
+ m_tree->Branch("err_ePHI_flt", &m_err_ePHI[1]);
+ m_tree->Branch("err_eTHETA_flt", &m_err_eTHETA[1]);
+ m_tree->Branch("err_eQOP_flt", &m_err_eQOP[1]);
+ m_tree->Branch("err_eT_flt", &m_err_eT[1]);
+ m_tree->Branch("pull_eLOC0_flt", &m_pull_eLOC0[1]);
+ m_tree->Branch("pull_eLOC1_flt", &m_pull_eLOC1[1]);
+ m_tree->Branch("pull_ePHI_flt", &m_pull_ePHI[1]);
+ m_tree->Branch("pull_eTHETA_flt", &m_pull_eTHETA[1]);
+ m_tree->Branch("pull_eQOP_flt", &m_pull_eQOP[1]);
+ m_tree->Branch("pull_eT_flt", &m_pull_eT[1]);
+ m_tree->Branch("g_x_flt", &m_x[1]);
+ m_tree->Branch("g_y_flt", &m_y[1]);
+ m_tree->Branch("g_z_flt", &m_z[1]);
+ m_tree->Branch("px_flt", &m_px[1]);
+ m_tree->Branch("py_flt", &m_py[1]);
+ m_tree->Branch("pz_flt", &m_pz[1]);
+ m_tree->Branch("eta_flt", &m_eta[1]);
+ m_tree->Branch("pT_flt", &m_pT[1]);
+
+ m_tree->Branch("nSmoothed", &m_nParams[2]);
+ m_tree->Branch("smoothed", &m_hasParams[2]);
+ m_tree->Branch("eLOC0_smt", &m_eLOC0[2]);
+ m_tree->Branch("eLOC1_smt", &m_eLOC1[2]);
+ m_tree->Branch("ePHI_smt", &m_ePHI[2]);
+ m_tree->Branch("eTHETA_smt", &m_eTHETA[2]);
+ m_tree->Branch("eQOP_smt", &m_eQOP[2]);
+ m_tree->Branch("eT_smt", &m_eT[2]);
+ m_tree->Branch("res_eLOC0_smt", &m_res_eLOC0[2]);
+ m_tree->Branch("res_eLOC1_smt", &m_res_eLOC1[2]);
+ m_tree->Branch("res_ePHI_smt", &m_res_ePHI[2]);
+ m_tree->Branch("res_eTHETA_smt", &m_res_eTHETA[2]);
+ m_tree->Branch("res_eQOP_smt", &m_res_eQOP[2]);
+ m_tree->Branch("res_eT_smt", &m_res_eT[2]);
+ m_tree->Branch("err_eLOC0_smt", &m_err_eLOC0[2]);
+ m_tree->Branch("err_eLOC1_smt", &m_err_eLOC1[2]);
+ m_tree->Branch("err_ePHI_smt", &m_err_ePHI[2]);
+ m_tree->Branch("err_eTHETA_smt", &m_err_eTHETA[2]);
+ m_tree->Branch("err_eQOP_smt", &m_err_eQOP[2]);
+ m_tree->Branch("err_eT_smt", &m_err_eT[2]);
+ m_tree->Branch("pull_eLOC0_smt", &m_pull_eLOC0[2]);
+ m_tree->Branch("pull_eLOC1_smt", &m_pull_eLOC1[2]);
+ m_tree->Branch("pull_ePHI_smt", &m_pull_ePHI[2]);
+ m_tree->Branch("pull_eTHETA_smt", &m_pull_eTHETA[2]);
+ m_tree->Branch("pull_eQOP_smt", &m_pull_eQOP[2]);
+ m_tree->Branch("pull_eT_smt", &m_pull_eT[2]);
+ m_tree->Branch("g_x_smt", &m_x[2]);
+ m_tree->Branch("g_y_smt", &m_y[2]);
+ m_tree->Branch("g_z_smt", &m_z[2]);
+ m_tree->Branch("px_smt", &m_px[2]);
+ m_tree->Branch("py_smt", &m_py[2]);
+ m_tree->Branch("pz_smt", &m_pz[2]);
+ m_tree->Branch("eta_smt", &m_eta[2]);
+ m_tree->Branch("pT_smt", &m_pT[2]);
+
+ m_tree->Branch("chi2", &m_chi2);
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode RootTrajectoryStatesWriterTool::finalize() {
+ m_file->cd();
+ m_tree->Write();
+ m_file->Close();
+ return StatusCode::SUCCESS;
+}
+
+StatusCode RootTrajectoryStatesWriterTool::write(
+ TrajectoriesContainer trajectories, Acts::GeometryContext gctx,
+ std::vector<Acts::CurvilinearTrackParameters> traj) const {
+
+ // 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);
+
+ // Get the event number
+ EventContext ctx = Gaudi::Hive::currentContext();
+ m_eventNr = ctx.eventID().event_number();
+
+ // Loop over the trajectories
+ for (size_t itraj = 0; itraj < trajectories.size(); ++itraj) {
+ const auto& traj = trajectories[itraj];
+
+ if (traj.empty()) {
+ ATH_MSG_WARNING("Empty trajectories object " << itraj);
+ continue;
+ }
+
+ // The trajectory index
+ m_multiTrajNr = itraj;
+
+ // The trajectory entry indices and the multiTrajectory
+ const auto& mj = traj.multiTrajectory();
+ const auto& trackTips = traj.tips();
+
+ // Loop over the entry indices for the subtrajectories
+ for (unsigned int isubtraj = 0; isubtraj < trackTips.size(); ++isubtraj) {
+ // The subtrajectory index
+ m_subTrajNr = isubtraj;
+ // The entry index for this subtrajectory
+ const auto& trackTip = trackTips[isubtraj];
+ // Collect the trajectory summary info
+ auto trajState =
+ Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
+ m_nMeasurements = trajState.nMeasurements;
+ m_nStates = trajState.nStates;
+
+ // Get the majority truth particle to this track
+ // int truthQ = 1.;
+ // identifyContributingParticles(hitParticlesMap, traj, trackTip,
+ // particleHitCounts);
+ // if (not particleHitCounts.empty()) {
+ // // Get the barcode of the majority truth particle
+ // auto barcode = particleHitCounts.front().particleId.value();
+ // // Find the truth particle via the barcode
+ // auto ip = particles.find(barcode);
+ // if (ip != particles.end()) {
+ // const auto& particle = *ip;
+ // ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
+ // // Get the truth particle charge
+ // truthQ = particle.charge();
+ // } else {
+ // ATH_MSG_WARNING("Truth particle with barcode = " << barcode
+ // << " not found!");
+ // }
+ // }
+
+ // Get the trackStates on the trajectory
+ m_nParams = {0, 0, 0};
+ mj.visitBackwards(trackTip, [&](const auto& state) {
+ // we only fill the track states with non-outlier measurement
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
+ return true;
+ }
+
+ const auto& surface = state.referenceSurface();
+
+ // get the truth hits corresponding to this trackState
+ // Use average truth in the case of multiple contributing sim hits
+ // const auto hitIdx = state.uncalibrated().index();
+ // auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
+ // auto [truthLocal, truthPos4, truthUnitDir] =
+ // averageSimHits(ctx.geoContext, surface, simHits, indices);
+ // // momemtum averaging makes even less sense than averaging position and
+ // // direction. use the first momentum or set q/p to zero
+ // float truthQOP = 0.0f;
+ // if (not indices.empty()) {
+ // // we assume that the indices are within valid ranges so we do not
+ // // need to check their validity again.
+ // const auto simHitIdx0 = indices.begin()->second;
+ // const auto& simHit0 = *simHits.nth(simHitIdx0);
+ // const auto p =
+ // simHit0.momentum4Before().template segment<3>(Acts::eMom0).norm();
+ // truthQOP = truthQ / p;
+ // }
+
+ // // fill the truth hit info
+ // m_t_x.push_back(truthPos4[Acts::ePos0]);
+ // m_t_y.push_back(truthPos4[Acts::ePos1]);
+ // m_t_z.push_back(truthPos4[Acts::ePos2]);
+ // m_t_r.push_back(perp(truthPos4.template segment<3>(Acts::ePos0)));
+ // 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 = truthPos4[Acts::eTime];
+ // 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();
+ m_volumeID.push_back(geoID.volume());
+ m_layerID.push_back(geoID.layer());
+ m_moduleID.push_back(geoID.sensitive());
+
+ // get the path length
+ m_pathLength.push_back(state.pathLength());
+
+ // expand the local measurements into the full bound space
+ Acts::BoundVector meas =
+ state.projector().transpose() * state.calibrated();
+ // extract local and global position
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ Acts::Vector3 mom(1, 1, 1);
+ Acts::Vector3 global =
+ surface.localToGlobal(gctx, local, mom);
+
+ // fill the measurement info
+ m_lx_hit.push_back(local[Acts::ePos0]);
+ m_ly_hit.push_back(local[Acts::ePos1]);
+ m_x_hit.push_back(global[Acts::ePos0]);
+ m_y_hit.push_back(global[Acts::ePos1]);
+ m_z_hit.push_back(global[Acts::ePos2]);
+
+ // status of the fitted track parameters
+ std::array<bool, 3> hasParams = {false, false, false};
+ // optional fitted track parameters
+ std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>>
+ trackParamsOpt = std::nullopt;
+ // lambda to get the fitted track parameters
+ auto getTrackParams = [&](unsigned int ipar) {
+ if (ipar == 0 && state.hasPredicted()) {
+ hasParams[0] = true;
+ m_nParams[0]++;
+ trackParamsOpt =
+ std::make_pair(state.predicted(), state.predictedCovariance());
+ } else if (ipar == 1 && state.hasFiltered()) {
+ hasParams[1] = true;
+ m_nParams[1]++;
+ trackParamsOpt =
+ std::make_pair(state.filtered(), state.filteredCovariance());
+ } else if (ipar == 2 && state.hasSmoothed()) {
+ hasParams[2] = true;
+ m_nParams[2]++;
+ trackParamsOpt =
+ std::make_pair(state.smoothed(), state.smoothedCovariance());
+ }
+ };
+
+ // fill the fitted track parameters
+ for (unsigned int ipar = 0; ipar < 3; ++ipar) {
+ // get the fitted track parameters
+ getTrackParams(ipar);
+ if (trackParamsOpt) {
+ const auto& [parameters, covariance] = *trackParamsOpt;
+ if (ipar == 0) {
+ //
+ // local hit residual info
+ auto H = state.effectiveProjector();
+ auto resCov = state.effectiveCalibratedCovariance() +
+ H * covariance * H.transpose();
+ 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_err_x_hit.push_back(
+ sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ m_err_y_hit.push_back(
+ sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+ 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_dim_hit.push_back(state.calibratedSize());
+ }
+
+ // track parameters
+ m_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0]);
+ m_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1]);
+ m_ePHI[ipar].push_back(parameters[Acts::eBoundPhi]);
+ m_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta]);
+ m_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP]);
+ 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)));
+
+ // further track parameter info
+ Acts::FreeVector freeParams =
+ Acts::detail::transformBoundToFreeParameters(surface, gctx, parameters);
+ m_x[ipar].push_back(freeParams[Acts::eFreePos0]);
+ m_y[ipar].push_back(freeParams[Acts::eFreePos1]);
+ m_z[ipar].push_back(freeParams[Acts::eFreePos2]);
+ auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
+ m_px[ipar].push_back(p * freeParams[Acts::eFreeDir0]);
+ m_py[ipar].push_back(p * freeParams[Acts::eFreeDir1]);
+ m_pz[ipar].push_back(p * freeParams[Acts::eFreeDir2]);
+ m_pT[ipar].push_back(p * std::hypot(freeParams[Acts::eFreeDir0], freeParams[Acts::eFreeDir1]));
+ m_eta[ipar].push_back(Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
+ } 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.);
+ }
+ // 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.);
+ }
+ // fill the track parameters status
+ m_hasParams[ipar].push_back(hasParams[ipar]);
+ }
+
+ // fill the chi2
+ m_chi2.push_back(state.chi2());
+
+ return true;
+ }); // all states
+
+ // fill the variables for one track to tree
+ m_tree->Fill();
+
+ // now reset
+ m_t_x.clear();
+ m_t_y.clear();
+ m_t_z.clear();
+ m_t_r.clear();
+ m_t_dx.clear();
+ m_t_dy.clear();
+ m_t_dz.clear();
+ m_t_eLOC0.clear();
+ m_t_eLOC1.clear();
+ m_t_ePHI.clear();
+ m_t_eTHETA.clear();
+ m_t_eQOP.clear();
+ m_t_eT.clear();
+
+ m_volumeID.clear();
+ m_layerID.clear();
+ m_moduleID.clear();
+ m_pathLength.clear();
+ m_lx_hit.clear();
+ m_ly_hit.clear();
+ m_x_hit.clear();
+ m_y_hit.clear();
+ m_z_hit.clear();
+ m_res_x_hit.clear();
+ m_res_y_hit.clear();
+ m_err_x_hit.clear();
+ m_err_y_hit.clear();
+ m_pull_x_hit.clear();
+ m_pull_y_hit.clear();
+ m_dim_hit.clear();
+
+ for (unsigned int ipar = 0; ipar < 3; ++ipar) {
+ m_hasParams[ipar].clear();
+ m_eLOC0[ipar].clear();
+ m_eLOC1[ipar].clear();
+ m_ePHI[ipar].clear();
+ m_eTHETA[ipar].clear();
+ m_eQOP[ipar].clear();
+ m_eT[ipar].clear();
+ m_res_eLOC0[ipar].clear();
+ m_res_eLOC1[ipar].clear();
+ m_res_ePHI[ipar].clear();
+ m_res_eTHETA[ipar].clear();
+ m_res_eQOP[ipar].clear();
+ m_res_eT[ipar].clear();
+ m_err_eLOC0[ipar].clear();
+ m_err_eLOC1[ipar].clear();
+ m_err_ePHI[ipar].clear();
+ m_err_eTHETA[ipar].clear();
+ m_err_eQOP[ipar].clear();
+ m_err_eT[ipar].clear();
+ m_pull_eLOC0[ipar].clear();
+ m_pull_eLOC1[ipar].clear();
+ m_pull_ePHI[ipar].clear();
+ m_pull_eTHETA[ipar].clear();
+ m_pull_eQOP[ipar].clear();
+ m_pull_eT[ipar].clear();
+ m_x[ipar].clear();
+ m_y[ipar].clear();
+ m_z[ipar].clear();
+ m_px[ipar].clear();
+ m_py[ipar].clear();
+ m_pz[ipar].clear();
+ m_eta[ipar].clear();
+ m_pT[ipar].clear();
+ }
+
+ m_chi2.clear();
+ } // all subtrajectories
+ } // all trajectories
+
+ return StatusCode::SUCCESS;
+}