Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Update athena dump reader #4011

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open

feat: Update athena dump reader #4011

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5405865
update athena dump reader
benjaminhuth Jan 9, 2025
12cd786
add comment
benjaminhuth Jan 9, 2025
8523ab1
fix include
benjaminhuth Jan 9, 2025
ac2b97e
update
benjaminhuth Jan 9, 2025
cee908b
update
benjaminhuth Jan 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 13 additions & 2 deletions Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ class RootAthenaDumpReader : public IReader {
// Name of tree
std::string treename;
// Name of inputfile
std::string inputfile;
std::vector<std::string> inputfiles;
// name of the output measurements
std::string outputMeasurements = "athena_measurements";
// name of the output pixel space points
Expand All @@ -58,14 +58,19 @@ class RootAthenaDumpReader : public IReader {
std::string outputClusters = "athena_clusters";
// name of the output particles
std::string outputParticles = "athena_particles";
// name of the simhit map
// name of the measurements -> particles map
std::string outputMeasurementParticlesMap = "athena_meas_parts_map";
// name of the particles -> measurements map
std::string outputParticleMeasurementsMap = "athena_parts_meas_map";
// name of the track parameters (fitted by athena?)
std::string outputTrackParameters = "athena_track_parameters";

/// Only extract spacepoints
bool onlySpacepoints = false;

/// Skip truth data
bool noTruth = false;

/// Only extract particles that passed the tracking requirements, for
/// details see:
/// https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetGNNTracking/src/DumpObjects.cxx?ref_type=heads#L1363
Expand Down Expand Up @@ -165,6 +170,8 @@ class RootAthenaDumpReader : public IReader {
this, "output_measurements"};
WriteDataHandle<IndexMultimap<ActsFatras::Barcode>> m_outputMeasParticleMap{
this, "output_meas_part_map"};
WriteDataHandle<InverseMultimap<ActsFatras::Barcode>> m_outputParticleMeasMap{
this, "output_part_meas_map"};

std::unique_ptr<const Acts::Logger> m_logger;
std::mutex m_read_mutex;
Expand Down Expand Up @@ -269,6 +276,9 @@ class RootAthenaDumpReader : public IReader {
std::vector<std::vector<float>> *SPstripCenterDistance{};
std::vector<std::vector<float>> *SPtopStripCenterPosition{};

// Those fields are not used currently
// Keep the code though, since it is annoying to write
/*
// Tracks
Int_t nTRK = 0;
Int_t TRKindex[maxTRK] = {}; //[nTRK]
Expand Down Expand Up @@ -299,5 +309,6 @@ class RootAthenaDumpReader : public IReader {
std::vector<std::vector<int>> *DTTstTruth_subDetType{};
std::vector<std::vector<int>> *DTTstTrack_subDetType{};
std::vector<std::vector<int>> *DTTstCommon_subDetType{};
*/
};
} // namespace ActsExamples
199 changes: 99 additions & 100 deletions Examples/Io/Root/src/RootAthenaDumpReader.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,8 @@ RootAthenaDumpReader::RootAthenaDumpReader(
: IReader(),
m_cfg(config),
m_logger(Acts::getDefaultLogger(name(), level)) {
if (m_cfg.inputfile.empty()) {
throw std::invalid_argument("Missing input filename");
if (m_cfg.inputfiles.empty()) {
throw std::invalid_argument("Empty input file list");
}
if (m_cfg.treename.empty()) {
throw std::invalid_argument("Missing tree name");
Expand All @@ -77,10 +77,13 @@ RootAthenaDumpReader::RootAthenaDumpReader(
m_outputStripSpacePoints.initialize(m_cfg.outputStripSpacePoints);
m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
if (!m_cfg.onlySpacepoints) {
m_outputClusters.initialize(m_cfg.outputClusters);
m_outputParticles.initialize(m_cfg.outputParticles);
m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap);
m_outputMeasurements.initialize(m_cfg.outputMeasurements);
if (!m_cfg.noTruth) {
m_outputClusters.initialize(m_cfg.outputClusters);
m_outputParticles.initialize(m_cfg.outputParticles);
m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap);
m_outputParticleMeasMap.initialize(m_cfg.outputParticleMeasurementsMap);
}
}

if (m_inputchain->GetBranch("SPtopStripDirection") == nullptr) {
Expand All @@ -89,35 +92,6 @@ RootAthenaDumpReader::RootAthenaDumpReader(
}

// Set the branches

// Set object pointer
CLhardware = nullptr;
CLparticleLink_eventIndex = nullptr;
CLparticleLink_barcode = nullptr;
CLbarcodesLinked = nullptr;
CLparticle_charge = nullptr;
CLphis = nullptr;
CLetas = nullptr;
CLtots = nullptr;
CLlocal_cov = nullptr;
Part_vParentID = nullptr;
Part_vParentBarcode = nullptr;
SPtopStripDirection = nullptr;
SPbottomStripDirection = nullptr;
SPstripCenterDistance = nullptr;
SPtopStripCenterPosition = nullptr;
TRKproperties = nullptr;
TRKpattern = nullptr;
TRKmeasurementsOnTrack_pixcl_sctcl_index = nullptr;
TRKoutliersOnTrack_pixcl_sctcl_index = nullptr;
TRKperigee_position = nullptr;
TRKperigee_momentum = nullptr;
DTTtrajectory_eventindex = nullptr;
DTTtrajectory_barcode = nullptr;
DTTstTruth_subDetType = nullptr;
DTTstTrack_subDetType = nullptr;
DTTstCommon_subDetType = nullptr;

m_inputchain->SetBranchAddress("run_number", &run_number);
m_inputchain->SetBranchAddress("event_number", &event_number);
m_inputchain->SetBranchAddress("nSE", &nSE);
Expand All @@ -134,12 +108,6 @@ RootAthenaDumpReader::RootAthenaDumpReader(
m_inputchain->SetBranchAddress("CLphi_module", CLphi_module);
m_inputchain->SetBranchAddress("CLside", CLside);
m_inputchain->SetBranchAddress("CLmoduleID", CLmoduleID);
m_inputchain->SetBranchAddress("CLparticleLink_eventIndex",
&CLparticleLink_eventIndex);
m_inputchain->SetBranchAddress("CLparticleLink_barcode",
&CLparticleLink_barcode);
m_inputchain->SetBranchAddress("CLbarcodesLinked", &CLbarcodesLinked);
m_inputchain->SetBranchAddress("CLparticle_charge", &CLparticle_charge);
m_inputchain->SetBranchAddress("CLphis", &CLphis);
m_inputchain->SetBranchAddress("CLetas", &CLetas);
m_inputchain->SetBranchAddress("CLtots", &CLtots);
Expand All @@ -161,28 +129,40 @@ RootAthenaDumpReader::RootAthenaDumpReader(
m_inputchain->SetBranchAddress("CLnorm_y", CLnorm_y);
m_inputchain->SetBranchAddress("CLnorm_z", CLnorm_z);
m_inputchain->SetBranchAddress("CLlocal_cov", &CLlocal_cov);
m_inputchain->SetBranchAddress("nPartEVT", &nPartEVT);
m_inputchain->SetBranchAddress("Part_event_number", Part_event_number);
m_inputchain->SetBranchAddress("Part_barcode", Part_barcode);
m_inputchain->SetBranchAddress("Part_px", Part_px);
m_inputchain->SetBranchAddress("Part_py", Part_py);
m_inputchain->SetBranchAddress("Part_pz", Part_pz);
m_inputchain->SetBranchAddress("Part_pt", Part_pt);
m_inputchain->SetBranchAddress("Part_eta", Part_eta);
m_inputchain->SetBranchAddress("Part_vx", Part_vx);
m_inputchain->SetBranchAddress("Part_vy", Part_vy);
m_inputchain->SetBranchAddress("Part_vz", Part_vz);
m_inputchain->SetBranchAddress("Part_radius", Part_radius);
m_inputchain->SetBranchAddress("Part_status", Part_status);
m_inputchain->SetBranchAddress("Part_charge", Part_charge);
m_inputchain->SetBranchAddress("Part_pdg_id", Part_pdg_id);
m_inputchain->SetBranchAddress("Part_passed", Part_passed);
m_inputchain->SetBranchAddress("Part_vProdNin", Part_vProdNin);
m_inputchain->SetBranchAddress("Part_vProdNout", Part_vProdNout);
m_inputchain->SetBranchAddress("Part_vProdStatus", Part_vProdStatus);
m_inputchain->SetBranchAddress("Part_vProdBarcode", Part_vProdBarcode);
m_inputchain->SetBranchAddress("Part_vParentID", &Part_vParentID);
m_inputchain->SetBranchAddress("Part_vParentBarcode", &Part_vParentBarcode);
if (!m_cfg.noTruth) {
m_inputchain->SetBranchAddress("CLparticleLink_eventIndex",
&CLparticleLink_eventIndex);
m_inputchain->SetBranchAddress("CLparticleLink_barcode",
&CLparticleLink_barcode);
m_inputchain->SetBranchAddress("CLbarcodesLinked", &CLbarcodesLinked);
m_inputchain->SetBranchAddress("CLparticle_charge", &CLparticle_charge);
}

if (!m_cfg.noTruth) {
Comment on lines +132 to +141
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is there a reason to separate these ifs?

m_inputchain->SetBranchAddress("nPartEVT", &nPartEVT);
m_inputchain->SetBranchAddress("Part_event_number", Part_event_number);
m_inputchain->SetBranchAddress("Part_barcode", Part_barcode);
m_inputchain->SetBranchAddress("Part_px", Part_px);
m_inputchain->SetBranchAddress("Part_py", Part_py);
m_inputchain->SetBranchAddress("Part_pz", Part_pz);
m_inputchain->SetBranchAddress("Part_pt", Part_pt);
m_inputchain->SetBranchAddress("Part_eta", Part_eta);
m_inputchain->SetBranchAddress("Part_vx", Part_vx);
m_inputchain->SetBranchAddress("Part_vy", Part_vy);
m_inputchain->SetBranchAddress("Part_vz", Part_vz);
m_inputchain->SetBranchAddress("Part_radius", Part_radius);
m_inputchain->SetBranchAddress("Part_status", Part_status);
m_inputchain->SetBranchAddress("Part_charge", Part_charge);
m_inputchain->SetBranchAddress("Part_pdg_id", Part_pdg_id);
m_inputchain->SetBranchAddress("Part_passed", Part_passed);
m_inputchain->SetBranchAddress("Part_vProdNin", Part_vProdNin);
m_inputchain->SetBranchAddress("Part_vProdNout", Part_vProdNout);
m_inputchain->SetBranchAddress("Part_vProdStatus", Part_vProdStatus);
m_inputchain->SetBranchAddress("Part_vProdBarcode", Part_vProdBarcode);
m_inputchain->SetBranchAddress("Part_vParentID", &Part_vParentID);
m_inputchain->SetBranchAddress("Part_vParentBarcode", &Part_vParentBarcode);
}

m_inputchain->SetBranchAddress("nSP", &nSP);
m_inputchain->SetBranchAddress("SPindex", SPindex);
m_inputchain->SetBranchAddress("SPx", SPx);
Expand All @@ -191,19 +171,25 @@ RootAthenaDumpReader::RootAthenaDumpReader(
m_inputchain->SetBranchAddress("SPCL1_index", SPCL1_index);
m_inputchain->SetBranchAddress("SPCL2_index", SPCL2_index);
m_inputchain->SetBranchAddress("SPisOverlap", SPisOverlap);
m_inputchain->SetBranchAddress("SPradius", SPradius);
m_inputchain->SetBranchAddress("SPcovr", SPcovr);
m_inputchain->SetBranchAddress("SPcovz", SPcovz);
m_inputchain->SetBranchAddress("SPhl_topstrip", SPhl_topstrip);
m_inputchain->SetBranchAddress("SPhl_botstrip", SPhl_botstrip);
m_inputchain->SetBranchAddress("SPtopStripDirection", &SPtopStripDirection);
m_inputchain->SetBranchAddress("SPbottomStripDirection",
&SPbottomStripDirection);
m_inputchain->SetBranchAddress("SPstripCenterDistance",
&SPstripCenterDistance);
m_inputchain->SetBranchAddress("SPtopStripCenterPosition",
&SPtopStripCenterPosition);

if (m_haveStripFeatures) {
m_inputchain->SetBranchAddress("SPradius", SPradius);
m_inputchain->SetBranchAddress("SPcovr", SPcovr);
m_inputchain->SetBranchAddress("SPcovz", SPcovz);
m_inputchain->SetBranchAddress("SPhl_topstrip", SPhl_topstrip);
m_inputchain->SetBranchAddress("SPhl_botstrip", SPhl_botstrip);
m_inputchain->SetBranchAddress("SPtopStripDirection", &SPtopStripDirection);
m_inputchain->SetBranchAddress("SPbottomStripDirection",
&SPbottomStripDirection);
m_inputchain->SetBranchAddress("SPstripCenterDistance",
&SPstripCenterDistance);
m_inputchain->SetBranchAddress("SPtopStripCenterPosition",
&SPtopStripCenterPosition);
}

// These quantities are not used currently and thus commented out
// I would like to keep the code, since it is always a pain to write it
/*
m_inputchain->SetBranchAddress("nTRK", &nTRK);
m_inputchain->SetBranchAddress("TRKindex", TRKindex);
m_inputchain->SetBranchAddress("TRKtrack_fitter", TRKtrack_fitter);
Expand Down Expand Up @@ -239,9 +225,12 @@ RootAthenaDumpReader::RootAthenaDumpReader(
&DTTstTrack_subDetType);
m_inputchain->SetBranchAddress("DTTstCommon_subDetType",
&DTTstCommon_subDetType);
*/

m_inputchain->Add(m_cfg.inputfile.c_str());
ACTS_DEBUG("Adding file " << m_cfg.inputfile << " to tree" << m_cfg.treename);
for (const auto& file : m_cfg.inputfiles) {
m_inputchain->Add(file.c_str());
ACTS_DEBUG("Adding file '" << file << "' to tree " << m_cfg.treename);
Comment on lines +231 to +232
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nitpick: in case of an error this might be better

Suggested change
m_inputchain->Add(file.c_str());
ACTS_DEBUG("Adding file '" << file << "' to tree " << m_cfg.treename);
ACTS_DEBUG("Adding file '" << file << "' to tree " << m_cfg.treename);
m_inputchain->Add(file.c_str());

}

m_events = m_inputchain->GetEntries();

Expand Down Expand Up @@ -447,25 +436,26 @@ RootAthenaDumpReader::readMeasurements(
}

std::size_t measIndex = measurements.size();
imIdxMap.emplace(im, measIndex);
createMeasurement(measurements, geoId, digiPars);

// Create measurement particles map and particles container
for (const auto& [subevt, barcode] :
Acts::zip(CLparticleLink_eventIndex->at(im),
CLparticleLink_barcode->at(im))) {
auto dummyBarcode = concatInts(barcode, subevt);
// If we don't find the particle, create one with default values
if (particles.find(dummyBarcode) == particles.end()) {
ACTS_VERBOSE("Particle with subevt " << subevt << ", barcode "
<< barcode
<< "not found, create dummy one");
particles.emplace(dummyBarcode, Acts::PdgParticle::eInvalid);
if (!m_cfg.noTruth) {
// Create measurement particles map and particles container
for (const auto& [subevt, barcode] :
Acts::zip(CLparticleLink_eventIndex->at(im),
CLparticleLink_barcode->at(im))) {
auto dummyBarcode = concatInts(barcode, subevt);
// If we don't find the particle, create one with default values
if (particles.find(dummyBarcode) == particles.end()) {
ACTS_VERBOSE("Particle with subevt "
<< subevt << ", barcode " << barcode
<< "not found, create dummy one");
particles.emplace(dummyBarcode, Acts::PdgParticle::eInvalid);
}
measPartMap.insert(
std::pair<Index, ActsFatras::Barcode>{measIndex, dummyBarcode});
}
measPartMap.insert(
std::pair<Index, ActsFatras::Barcode>{measIndex, dummyBarcode});
}

imIdxMap.emplace(im, measIndex);
}

if (measurements.size() < static_cast<std::size_t>(nCL)) {
Expand All @@ -479,8 +469,8 @@ RootAthenaDumpReader::readMeasurements(
}

if (nTotalTotZero > 0) {
ACTS_WARNING(nTotalTotZero << " / " << nCL
<< " clusters have zero time-over-threshold");
ACTS_DEBUG(nTotalTotZero << " / " << nCL
<< " clusters have zero time-over-threshold");
}

return {clusters, measurements, measPartMap, imIdxMap};
Expand Down Expand Up @@ -618,7 +608,8 @@ RootAthenaDumpReader::readSpacepoints(
ACTS_DEBUG("Skipped " << skippedSpacePoints
<< " because of eta/phi overlaps");
}
if (spacePoints.size() < static_cast<std::size_t>(nSP)) {
if (spacePoints.size() <
(static_cast<std::size_t>(nSP) - skippedSpacePoints)) {
ACTS_WARNING("Could not convert " << nSP - spacePoints.size() << " of "
<< nSP << " spacepoints");
}
Expand Down Expand Up @@ -701,19 +692,27 @@ ProcessCode RootAthenaDumpReader::read(const AlgorithmContext& ctx) {
std::optional<std::unordered_map<int, std::size_t>> optImIdxMap;

if (!m_cfg.onlySpacepoints) {
auto candidateParticles = readParticles();
SimParticleContainer candidateParticles;

if (!m_cfg.noTruth) {
candidateParticles = readParticles();
}

auto [clusters, measurements, candidateMeasPartMap, imIdxMap] =
readMeasurements(candidateParticles, ctx.geoContext);
optImIdxMap.emplace(std::move(imIdxMap));

auto [particles, measPartMap] =
reprocessParticles(candidateParticles, candidateMeasPartMap);

m_outputClusters(ctx, std::move(clusters));
m_outputParticles(ctx, std::move(particles));
m_outputMeasParticleMap(ctx, std::move(measPartMap));
m_outputMeasurements(ctx, std::move(measurements));

if (!m_cfg.noTruth) {
auto [particles, measPartMap] =
reprocessParticles(candidateParticles, candidateMeasPartMap);

m_outputParticles(ctx, std::move(particles));
m_outputParticleMeasMap(ctx, invertIndexMultimap(measPartMap));
m_outputMeasParticleMap(ctx, std::move(measPartMap));
}
}

auto [spacePoints, pixelSpacePoints, stripSpacePoints] =
Expand Down
9 changes: 5 additions & 4 deletions Examples/Python/src/Input.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,11 +92,12 @@ void addInput(Context& ctx) {

ACTS_PYTHON_DECLARE_READER(
ActsExamples::RootAthenaDumpReader, mex, "RootAthenaDumpReader", treename,
inputfile, outputMeasurements, outputPixelSpacePoints,
inputfiles, outputMeasurements, outputPixelSpacePoints,
outputStripSpacePoints, outputSpacePoints, outputClusters,
outputMeasurementParticlesMap, outputParticles, onlyPassedParticles,
skipOverlapSPsPhi, skipOverlapSPsEta, geometryIdMap, trackingGeometry,
absBoundaryTolerance);
outputMeasurementParticlesMap, outputParticleMeasurementsMap,
outputParticles, onlyPassedParticles, skipOverlapSPsPhi,
skipOverlapSPsEta, geometryIdMap, trackingGeometry, absBoundaryTolerance,
onlySpacepoints, noTruth);

#ifdef WITH_GEOMODEL_PLUGIN
ACTS_PYTHON_DECLARE_READER(ActsExamples::RootAthenaDumpGeoIdCollector, mex,
Expand Down
Loading