From 71722df9b7e0af490d79a479173f57b7f59c7417 Mon Sep 17 00:00:00 2001
From: Xiaocong Ai <xiaocong.ai@cern.ch>
Date: Tue, 29 Oct 2024 09:57:53 +0100
Subject: [PATCH] fix truth barcode retrieving
---
.../src/RootTrajectoryStatesWriterTool.cxx | 11 ++++---
.../src/RootTrajectorySummaryWriterTool.cxx | 33 ++++++++++++++++++-
.../src/RootTrajectorySummaryWriterTool.h | 3 ++
.../src/TrackClassification.cxx | 12 +++----
4 files changed, 48 insertions(+), 11 deletions(-)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
index 785089cdb..a4f86a5a7 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
@@ -279,7 +279,6 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
const EventContext& ctx = Gaudi::Hive::currentContext();
m_eventNr = ctx.eventID().event_number();
- std::vector<ParticleHitCount> particleHitCounts;
std::map<int, const HepMC::GenParticle*> particles {};
std::map<std::pair<int, Identifier>, const FaserSiHit*> siHitMap;
@@ -294,7 +293,7 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
}
ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
for (const auto& particle : mcEvents->front()->particles()) {
- particles[particle->id()] = &(*particle);
+ particles[HepMC::barcode(particle)] = &(*particle);
}
SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
@@ -304,12 +303,13 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
ATH_CHECK(siHitCollection.isValid());
for (const FaserSiHit& hit : *siHitCollection) {
- int barcode = hit.trackNumber();
+ int trackNumber = 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;
+ siHitMap[std::make_pair(trackNumber, id)] = &hit;
}
}
+ std::vector<ParticleHitCount> particleHitCounts = {};
for (const auto& track : tracks) {
m_trackNr = track.index();
@@ -425,6 +425,9 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
m_t_eQOP.push_back(truthQOP);
m_t_eT.push_back(truthTIME);
} else {
+ if(isMC){
+ ATH_MSG_WARNING("Hit with particle barcode " << barcode <<" and wafer Id " << waferId <<" not found.");
+ }
m_t_x.push_back(nan);
m_t_y.push_back(nan);
m_t_z.push_back(nan);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
index 4ef9c4dda..e6323c8de 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.cxx
@@ -68,12 +68,15 @@ StatusCode RootTrajectorySummaryWriterTool::initialize() {
m_outputTree->Branch("nSharedHits", &m_nSharedHits);
m_outputTree->Branch("chi2Sum", &m_chi2Sum);
m_outputTree->Branch("NDF", &m_NDF);
+ m_outputTree->Branch("rMax_hit", &m_rMax_hit);
m_outputTree->Branch("measurementChi2", &m_measurementChi2);
m_outputTree->Branch("outlierChi2", &m_outlierChi2);
m_outputTree->Branch("measurementVolume", &m_measurementVolume);
m_outputTree->Branch("measurementLayer", &m_measurementLayer);
+ m_outputTree->Branch("measurementIndex", &m_measurementIndex);
m_outputTree->Branch("outlierVolume", &m_outlierVolume);
m_outputTree->Branch("outlierLayer", &m_outlierLayer);
+ m_outputTree->Branch("outlierIndex", &m_outlierIndex);
m_outputTree->Branch("nMajorityHits", &m_nMajorityHits);
m_outputTree->Branch("majorityParticleId", &m_majorityParticleId);
@@ -158,7 +161,7 @@ StatusCode RootTrajectorySummaryWriterTool::write(
}
for (const auto& particle : mcEvents->front()->particles()) {
- particles[particle->id()] = &(*particle);
+ particles[HepMC::barcode(particle)] = &(*particle);
}
}
@@ -168,6 +171,7 @@ StatusCode RootTrajectorySummaryWriterTool::write(
// Get the event number
m_eventNr = ctx.eventID().event_number();
+ Acts::Vector3 truthUnitDir(1,1,1);
// Loop over the trajectories
for (const auto& track : tracks) {
m_trackNr.push_back(track.index());
@@ -185,9 +189,12 @@ StatusCode RootTrajectorySummaryWriterTool::write(
std::vector<double> measurementChi2;
std::vector<double> measurementVolume;
std::vector<double> measurementLayer;
+ std::vector<int> measurementIndex;
std::vector<double> outlierChi2;
std::vector<double> outlierVolume;
std::vector<double> outlierLayer;
+ std::vector<int> outlierIndex;
+ double rMax_hit = 0;
for (const auto& state : track.trackStatesReversed()) {
const auto& geoID = state.referenceSurface().geometryId();
const auto& volume = geoID.volume();
@@ -196,22 +203,43 @@ StatusCode RootTrajectorySummaryWriterTool::write(
measurementChi2.push_back(state.chi2());
measurementVolume.push_back(volume);
measurementLayer.push_back(layer);
+ // expand the local measurements into the full bound space
+ Acts::BoundVector meas = state.effectiveProjector().transpose() *
+ state.effectiveCalibrated();
+ // extract local and global position
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ const auto& surface = state.referenceSurface();
+ // The direction does not matter for plane surface
+ Acts::Vector3 global =
+ surface.localToGlobal(geoContext, local, truthUnitDir);
+ auto r_hit = std::hypot(global[Acts::ePos0],global[Acts::ePos1]);
+ if(r_hit>rMax_hit){
+ rMax_hit = r_hit;
+ }
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ measurementIndex.push_back(sl.index());
}
if (state.typeFlags().test(Acts::TrackStateFlag::OutlierFlag)) {
outlierChi2.push_back(state.chi2());
outlierVolume.push_back(volume);
outlierLayer.push_back(layer);
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ outlierIndex.push_back(sl.index());
}
}
+
// IDs are stored as double (as the vector of vector of int is not known
// to ROOT)
m_measurementChi2.push_back(std::move(measurementChi2));
m_measurementVolume.push_back(std::move(measurementVolume));
m_measurementLayer.push_back(std::move(measurementLayer));
+ m_measurementIndex.push_back(std::move(measurementIndex));
m_outlierChi2.push_back(std::move(outlierChi2));
m_outlierVolume.push_back(std::move(outlierVolume));
m_outlierLayer.push_back(std::move(outlierLayer));
+ m_outlierIndex.push_back(std::move(outlierIndex));
+ m_rMax_hit.push_back(rMax_hit);
}
// Initialize the truth particle info
@@ -422,8 +450,11 @@ StatusCode RootTrajectorySummaryWriterTool::write(
m_outlierChi2.clear();
m_measurementVolume.clear();
m_measurementLayer.clear();
+ m_measurementIndex.clear();
m_outlierVolume.clear();
m_outlierLayer.clear();
+ m_outlierIndex.clear();
+ m_rMax_hit.clear();
m_nMajorityHits.clear();
m_majorityParticleId.clear();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.h
index 383431d78..91dcfbd53 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectorySummaryWriterTool.h
@@ -58,13 +58,16 @@ private:
mutable std::vector<unsigned int> m_nSharedHits; ///< The number of shared hits
mutable std::vector<float> m_chi2Sum; ///< The total chi2
mutable std::vector<unsigned int> m_NDF; ///< The number of ndf of the measurements+outliers
+ mutable std::vector<double> m_rMax_hit; ///< max radius of the measurements on the track
mutable std::vector<std::vector<double>> m_measurementChi2; ///< The chi2 on all measurement states
mutable std::vector<std::vector<double>> m_outlierChi2; ///< The chi2 on all outlier states
// FIXME replace volume, layer, ... with station, layer, ...
mutable std::vector<std::vector<double>> m_measurementVolume; ///< The volume id of the measurements
mutable std::vector<std::vector<double>> m_measurementLayer; ///< The layer id of the measurements
+ mutable std::vector<std::vector<int>> m_measurementIndex; ///< The hit id of the measurements
mutable std::vector<std::vector<double>> m_outlierVolume; ///< The volume id of the outliers
mutable std::vector<std::vector<double>> m_outlierLayer; ///< The layer id of the outliers
+ mutable std::vector<std::vector<int>> m_outlierIndex; ///< The hit id of the outliers
// The majority truth particle info
mutable std::vector<unsigned int> m_nMajorityHits; ///< The number of hits from majority particle
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
index a3cd8c500..ce2454ee5 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackClassification.cxx
@@ -52,17 +52,17 @@ void identifyContributingParticles(
IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
if(sl.hit()==nullptr){
throw std::runtime_error("The source hit is empty");
- }
+ }
for (const Identifier &id : sl.hit()->rdoList()) {
if (simDataCollection.count(id) == 0) {
- continue; //@todo:is this correct?
+ continue;
}
const auto &deposits = simDataCollection.find(id)->second.getdeposits();
for (const TrackerSimData::Deposit &deposit : deposits) {
- int barcode = deposit.first->id();
+ int barcode = deposit.first.barcode();
if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
barcodes.push_back(barcode);
- increaseHitCount(particleHitCounts, deposit.first->id());
+ increaseHitCount(particleHitCounts, barcode);
}
}
}
@@ -85,10 +85,10 @@ void identifyContributingParticles(
if (simDataCollection.count(id) == 0) continue;
const auto &deposits = simDataCollection.find(id)->second.getdeposits();
for (const TrackerSimData::Deposit &deposit : deposits) {
- int barcode = deposit.first->id();
+ int barcode = deposit.first.barcode();
if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
barcodes.push_back(barcode);
- increaseHitCount(particleHitCounts, deposit.first->id());
+ increaseHitCount(particleHitCounts, barcode);
}
}
}
--
GitLab