From d4c9de56ef8e09e61b91961555e955bb39df85a6 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 4 Sep 2026 08:15:16 +0200 Subject: [PATCH 1/3] Make the geometry doctor's reachability audit trustworthy This fixes two defects in the reachability audit, adds a second measurement to it, and raises its default sampling. - TGeoManager::FindNode() resumes from the navigator's current branch, so a walk that samples one placement at a time asked every question from inside the placement it was testing, and that placement won wherever two volumes overlap. A control geometry with two mutually overlapping boxes had both of them reporting themselves 100% reached. Each query now starts from CdTop(). - The path test was a bare rfind at offset 0, so .../X_1 matched .../X_10 and a point the navigator gave to a different sibling counted as reached. 155 of the 1570 mothers in the ALICE geometry have sibling names where one is a prefix of another, barrel among them with SMOD_1 and SMOD_10. - Reachability counts a point that lands in a daughter, which is the wrong question for material. For every point that is nominally the volume's own medium, inside its shape and inside none of its daughters, FindNode must now return exactly that path; the share that does is reported as "own kept". The daughter test is the insideAnyDaughter() the field classification already uses. - The default sample count goes from 32 to 1000, and the list is sorted worst first. At 32 a single boundary point crosses the 0.999 threshold. On the full ALICE geometry at 2000 samples the HMPID absorbers move from 99.3% reached to 5.5% and 8.6%, the FT0 mirror strips to 24%, the ITS cage foam to 36%; 391 placements are partially shadowed rather than 3, and one more is never reached at all, CageEndCap_1/CageEndCapRoundCross_1, which ROOT's own CheckOverlaps also reports. The corrected absorber figures agree with a direct containment measurement, in which 94% and 90% of the two plates lie inside the B077 space-frame envelope. The audit costs 4.7x more for this, 80 s instead of 17 s at 1000 samples on one core. Co-Authored-By: Claude Opus 5 --- run/o2sim_geometry_doctor.cxx | 108 ++++++++++++++++++++++++++-------- 1 file changed, 84 insertions(+), 24 deletions(-) diff --git a/run/o2sim_geometry_doctor.cxx b/run/o2sim_geometry_doctor.cxx index 1664580766276..2904aaf6d4bb3 100644 --- a/run/o2sim_geometry_doctor.cxx +++ b/run/o2sim_geometry_doctor.cxx @@ -718,14 +718,41 @@ bool supportFromJson(const json& in, Support& support) /// common defect lives -- a daughter outside its mother is a property of the node, /// not of the path that reaches it. A node whose mother is itself placed many /// times is therefore sampled once, in the first of those placements. +/// +/// Two questions are asked of every point, and they are not the same question. +/// *Reachability* asks whether the navigator's path passes through this placement +/// at all, so a point that lands in one of its own daughters counts. *Self +/// material* asks the stronger question: for a point that is nominally this +/// volume's own material -- inside its shape and inside none of its daughters -- +/// FindNode() must return exactly this path, not a prefix of it and not something +/// else. A mother whose own medium is entirely taken by an overlapping foreign +/// volume is still "reached" through its daughters, and only the second question +/// sees that its material is gone. struct Reach { std::string medium, mother, worstPath; long sampled = 0; double fraction = 1.; + long ownSampled = 0; ///< points that are nominally this volume's own material + double ownFraction = 1.; ///< of those, the share the navigator actually gives it }; constexpr int kReachRejectionTries = 400; +/// `found` passes through `path` -- a prefix match that must end on a path +/// separator. Without the boundary check `.../X_1` matches `.../X_10`, and copy +/// numbers 1 and 10 in one mother are common enough in ALICE that the check would +/// silently accept a point the navigator gave to a different sibling. +inline bool passesThrough(const std::string& found, const std::string& path) +{ + return found.compare(0, path.size(), path) == 0 && + (found.size() == path.size() || found[path.size()] == '/'); +} + +/// Defined with the placement table below. Deliberately the same predicate the +/// field classification already uses for "own material", so the two parts of this +/// tool cannot disagree about what a volume's own material is. +bool insideAnyDaughter(TGeoVolume* volume, const double* local); + class ReachAudit { public: @@ -785,37 +812,61 @@ void ReachAudit::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::stri // an assembly is expanded away at closure, so FindNode never returns one if (!volume->IsAssembly()) { - int drawn = 0, reached = 0; + int drawn = 0, reached = 0, ownDrawn = 0, ownReached = 0; + const bool hasDaughters = volume->GetNdaughters() > 0; for (int i = 0; i < mSamples; ++i) { double local[3], global[3]; if (!samplePoint(volume->GetShape(), local)) { break; } ++drawn; + // nominally this volume's own material: inside its shape, inside none of its + // daughters. A leaf owns every point of its shape, so skip the walk there. + const bool own = !hasDaughters || !insideAnyDaughter(volume, local); + if (own) { + ++ownDrawn; + } here.LocalToMaster(local, global); + // FindNode() resumes from wherever the navigator currently is, so without + // this the audit asks each question from inside the very placement it is + // testing and that placement wins every genuinely ambiguous point. Two + // mutually overlapping volumes then both report themselves fully reached. + // Starting from the top makes the answer the navigator's own, and the same + // one a track crossing the region would get. + gGeoManager->CdTop(); if (gGeoManager->FindNode(global[0], global[1], global[2]) == nullptr) { continue; } const std::string found = gGeoManager->GetPath(); - // reached if the navigator's own path passes through this placement - if (found.rfind(myPath, 0) == 0) { - ++reached; + if (!passesThrough(found, myPath)) { + continue; + } + ++reached; // the navigator's own path passes through this placement + if (own && found.size() == myPath.size()) { + ++ownReached; // ...and it stopped here, so the material really is this one's } } if (drawn == 0) { ++mUnsampleable; // a sliver too thin for the rejection budget; says nothing } else { ++mSampled; - const double fraction = double(reached) / drawn; - if (fraction < 0.999) { - auto* medium = volume->GetMedium(); - Reach entry; - entry.medium = medium != nullptr ? medium->GetName() : "(none)"; - entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; - entry.worstPath = myPath; - entry.sampled = drawn; - entry.fraction = fraction; - (fraction == 0. ? mDead : mPartial).push_back(entry); + auto* medium = volume->GetMedium(); + Reach entry; + entry.medium = medium != nullptr ? medium->GetName() : "(none)"; + entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; + entry.worstPath = myPath; + entry.sampled = drawn; + entry.fraction = double(reached) / drawn; + entry.ownSampled = ownDrawn; + entry.ownFraction = ownDrawn > 0 ? double(ownReached) / ownDrawn : 1.; + // Either number can fail on its own. A mother almost entirely filled by its + // daughters keeps a high reached fraction while the sliver of its own medium + // is taken by a foreign volume, and that sliver is the material that + // disappears -- so classify on whichever of the two is worse. + if (entry.fraction == 0.) { + mDead.push_back(entry); + } else if (entry.fraction < 0.999 || entry.ownFraction < 0.999) { + mPartial.push_back(entry); } } } @@ -839,6 +890,11 @@ long reportReachability(int samples, Report& report) audit.nodesVisited(), audit.nodesSampled(), audit.nodesUnsampleable())); report(form(" %ld placements the navigator never reaches, %zu it reaches only in part", (long)audit.dead().size(), audit.partial().size())); + // worst first: with hundreds of small overlaps the walk order is not a ranking + auto partial = audit.partial(); + std::sort(partial.begin(), partial.end(), [](const Reach& a, const Reach& b) { + return std::min(a.fraction, a.ownFraction) < std::min(b.fraction, b.ownFraction); + }); if (!audit.dead().empty()) { report(" unreachable -- these carry no material and produce no hits:"); report(form(" %-12s %-18s %10s %s", "mother", "medium", "sampled", "path")); @@ -847,17 +903,20 @@ long reportReachability(int samples, Report& report) entry.worstPath.c_str())); } } - for (size_t i = 0; i < audit.partial().size() && i < 20; ++i) { - const auto& entry = audit.partial()[i]; + for (size_t i = 0; i < partial.size() && i < 20; ++i) { + const auto& entry = partial[i]; if (i == 0) { - report(" partially shadowed -- an overlapping sibling or an extruding placement:"); - report(form(" %-12s %-18s %8s %s", "mother", "medium", "reached", "path")); + report(" partially shadowed -- an overlapping sibling or an extruding placement."); + report(" 'reached' is how much of the placement the navigator enters at all; 'own kept'"); + report(" how much of the medium this volume was given to carry survives as its own:"); + report(form(" %-12s %-18s %8s %9s %s", "mother", "medium", "reached", "own kept", "path")); } - report(form(" %-12s %-18s %7.1f%% %s", entry.mother.c_str(), entry.medium.c_str(), - 100. * entry.fraction, entry.worstPath.c_str())); + report(form(" %-12s %-18s %7.1f%% %8.1f%% %s", entry.mother.c_str(), entry.medium.c_str(), + 100. * entry.fraction, 100. * entry.ownFraction, entry.worstPath.c_str())); } - if (audit.partial().size() > 20) { - report(form(" ... and %zu more", audit.partial().size() - 20)); + if (partial.size() > 20) { + report(form(" ... and %zu more, all above %.1f%%", partial.size() - 20, + 100. * std::min(partial[19].fraction, partial[19].ownFraction))); } report(""); return (long)audit.dead().size(); @@ -1710,8 +1769,9 @@ int main(int argc, char** argv) "prefix for the report, the proposals and the placement table") // ("verify-anchors", bpo::value(&options.anchorFile), // "check the classification against known-good volumes listed in this JSON file") // - ("reachability-samples", bpo::value(&options.reachSamples)->default_value(32), // - "points drawn inside each placement for the reachability audit; 0 disables it") // + ("reachability-samples", bpo::value(&options.reachSamples)->default_value(1000), // + "points drawn inside each placement for the reachability audit; 0 disables it. Below a few " // + "hundred the audit reports genuine placements as partially shadowed") // ("reachability-only", bpo::bool_switch(&options.reachabilityOnly), // "run only the reachability audit, which needs no magnetic field"); From f9b28fcaf8bee73934197655f4388c9491f7e9d8 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 4 Sep 2026 09:21:05 +0200 Subject: [PATCH 2/3] Multi-thread the geometry doctor's reachability audit This makes the reachability audit use every core and makes its answer independent of how many it used. - The tree walk is separated from the sampling. The walk collects one representative placement per node object; the placements are then sampled in parallel, each thread with its own ROOT navigator. - Each placement draws from its own generator, seeded from its index. A single shared generator made every placement's numbers depend on the order the others were sampled in. - The navigator's answer is compared node by node rather than as a path string, which also drops the path formatting from the inner loop. - 1000 samples on the full ALICE geometry: 158 s on one core before, 68 s on one core now, 3.1 s on 128. The report is byte-identical at every thread count. Co-Authored-By: Claude Opus 5 --- run/o2sim_geometry_doctor.cxx | 319 ++++++++++++++++++++++------------ 1 file changed, 207 insertions(+), 112 deletions(-) diff --git a/run/o2sim_geometry_doctor.cxx b/run/o2sim_geometry_doctor.cxx index 2904aaf6d4bb3..87d41a4f83167 100644 --- a/run/o2sim_geometry_doctor.cxx +++ b/run/o2sim_geometry_doctor.cxx @@ -50,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -63,6 +64,7 @@ #include #include +#include #include #include #include @@ -72,6 +74,7 @@ #include #include #include +#include #include #include @@ -732,56 +735,81 @@ struct Reach { std::string medium, mother, worstPath; long sampled = 0; double fraction = 1.; - long ownSampled = 0; ///< points that are nominally this volume's own material + long ownSampled = 0; ///< points that are nominally this volume's own material double ownFraction = 1.; ///< of those, the share the navigator actually gives it }; constexpr int kReachRejectionTries = 400; -/// `found` passes through `path` -- a prefix match that must end on a path -/// separator. Without the boundary check `.../X_1` matches `.../X_10`, and copy -/// numbers 1 and 10 in one mother are common enough in ALICE that the check would -/// silently accept a point the navigator gave to a different sibling. -inline bool passesThrough(const std::string& found, const std::string& path) -{ - return found.compare(0, path.size(), path) == 0 && - (found.size() == path.size() || found[path.size()] == '/'); -} - /// Defined with the placement table below. Deliberately the same predicate the /// field classification already uses for "own material", so the two parts of this /// tool cannot disagree about what a volume's own material is. bool insideAnyDaughter(TGeoVolume* volume, const double* local); -class ReachAudit -{ - public: - explicit ReachAudit(int samples) : mSamples(samples) {} +/// One placement to sample: the node, where it sits, and the chain of nodes that +/// reaches it. The chain is what the navigator's answer is compared against -- +/// node identity rather than a path string, so no formatting or copy-number +/// ambiguity can enter the comparison. +struct ReachTask { + TGeoNode* node = nullptr; + TGeoHMatrix matrix; + std::string path; + std::vector chain; ///< top node first, this node last +}; - void walk(TGeoNode* node) { walk(node, TGeoHMatrix(), ""); } +struct ReachResult { + bool sampled = false; + long drawn = 0, reached = 0, ownDrawn = 0, ownReached = 0; +}; - const std::vector& dead() const { return mDead; } - const std::vector& partial() const { return mPartial; } +/// Collects one representative placement per node object. This half of the audit +/// is inherently sequential -- it carries the matrix chain down the tree and prunes +/// on node identity -- but it is also cheap, because it draws no points. +class ReachCollector +{ + public: + void walk(TGeoNode* node) { walk(node, TGeoHMatrix(), "", {}); } + std::vector& tasks() { return mTasks; } long nodesVisited() const { return mVisited; } - long nodesSampled() const { return mSampled; } - long nodesUnsampleable() const { return mUnsampleable; } private: - void walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path); - bool samplePoint(TGeoShape* shape, double* local); - - int mSamples; - long mVisited = 0, mSampled = 0, mUnsampleable = 0; - TRandom3 mRandom{20260901}; + void walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path, std::vector chain); std::set mSeen; - std::vector mDead, mPartial; + std::vector mTasks; + long mVisited = 0; }; +void ReachCollector::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path, + std::vector chain) +{ + if (!mSeen.insert(node).second) { + return; // this node object, and therefore its whole subtree, is already covered + } + TGeoHMatrix here = parent; + here.Multiply(node->GetMatrix()); + const std::string myPath = path + "/" + node->GetName(); + chain.push_back(node); + ++mVisited; + + // an assembly is expanded away at closure, so FindNode never returns one + if (!node->GetVolume()->IsAssembly()) { + ReachTask task; + task.node = node; + task.matrix = here; + task.path = myPath; + task.chain = chain; + mTasks.push_back(std::move(task)); + } + for (int i = 0; i < node->GetNdaughters(); ++i) { + walk(node->GetDaughter(i), here, myPath, chain); + } +} + /// Rejection sampling against the shape itself. A TGeoCompositeShape inherits /// TGeoBBox, so its DX/DY/DZ describe a box that still contains the holes and /// subtractions -- only Contains() knows the difference. GetOrigin() matters too: /// the box need not be centred on the local origin. -bool ReachAudit::samplePoint(TGeoShape* shape, double* local) +bool samplePoint(TGeoShape* shape, TRandom3& random, double* local) { auto* box = dynamic_cast(shape); if (box == nullptr) { @@ -789,9 +817,9 @@ bool ReachAudit::samplePoint(TGeoShape* shape, double* local) } const double* origin = box->GetOrigin(); for (int attempt = 0; attempt < kReachRejectionTries; ++attempt) { - local[0] = origin[0] + box->GetDX() * (2. * mRandom.Rndm() - 1.); - local[1] = origin[1] + box->GetDY() * (2. * mRandom.Rndm() - 1.); - local[2] = origin[2] + box->GetDZ() * (2. * mRandom.Rndm() - 1.); + local[0] = origin[0] + box->GetDX() * (2. * random.Rndm() - 1.); + local[1] = origin[1] + box->GetDY() * (2. * random.Rndm() - 1.); + local[2] = origin[2] + box->GetDZ() * (2. * random.Rndm() - 1.); if (shape->Contains(local)) { return true; } @@ -799,106 +827,169 @@ bool ReachAudit::samplePoint(TGeoShape* shape, double* local) return false; } -void ReachAudit::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path) +/// One generator per placement, seeded from its index. A single shared generator +/// would make every placement's numbers depend on the order the others were +/// sampled in -- which is the walk order in a serial run and nothing at all in a +/// parallel one. Seeding per placement makes the audit reproducible and identical +/// whatever --jobs is set to. +unsigned int seedFor(size_t index) { - if (!mSeen.insert(node).second) { - return; // this node object, and therefore its whole subtree, is already covered - } - TGeoHMatrix here = parent; - here.Multiply(node->GetMatrix()); - TGeoVolume* volume = node->GetVolume(); - const std::string myPath = path + "/" + node->GetName(); - ++mVisited; + unsigned long long x = 20260901ull + 0x9E3779B97F4A7C15ull * (index + 1); + x ^= x >> 30; + x *= 0xBF58476D1CE4E5B9ull; + x ^= x >> 27; + return (unsigned int)(x >> 33) | 1u; +} - // an assembly is expanded away at closure, so FindNode never returns one - if (!volume->IsAssembly()) { - int drawn = 0, reached = 0, ownDrawn = 0, ownReached = 0; - const bool hasDaughters = volume->GetNdaughters() > 0; - for (int i = 0; i < mSamples; ++i) { - double local[3], global[3]; - if (!samplePoint(volume->GetShape(), local)) { - break; - } - ++drawn; - // nominally this volume's own material: inside its shape, inside none of its - // daughters. A leaf owns every point of its shape, so skip the walk there. - const bool own = !hasDaughters || !insideAnyDaughter(volume, local); - if (own) { - ++ownDrawn; - } - here.LocalToMaster(local, global); - // FindNode() resumes from wherever the navigator currently is, so without - // this the audit asks each question from inside the very placement it is - // testing and that placement wins every genuinely ambiguous point. Two - // mutually overlapping volumes then both report themselves fully reached. - // Starting from the top makes the answer the navigator's own, and the same - // one a track crossing the region would get. - gGeoManager->CdTop(); - if (gGeoManager->FindNode(global[0], global[1], global[2]) == nullptr) { - continue; - } - const std::string found = gGeoManager->GetPath(); - if (!passesThrough(found, myPath)) { - continue; - } - ++reached; // the navigator's own path passes through this placement - if (own && found.size() == myPath.size()) { - ++ownReached; // ...and it stopped here, so the material really is this one's - } - } - if (drawn == 0) { - ++mUnsampleable; // a sliver too thin for the rejection budget; says nothing - } else { - ++mSampled; - auto* medium = volume->GetMedium(); - Reach entry; - entry.medium = medium != nullptr ? medium->GetName() : "(none)"; - entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; - entry.worstPath = myPath; - entry.sampled = drawn; - entry.fraction = double(reached) / drawn; - entry.ownSampled = ownDrawn; - entry.ownFraction = ownDrawn > 0 ? double(ownReached) / ownDrawn : 1.; - // Either number can fail on its own. A mother almost entirely filled by its - // daughters keeps a high reached fraction while the sliver of its own medium - // is taken by a foreign volume, and that sliver is the material that - // disappears -- so classify on whichever of the two is worse. - if (entry.fraction == 0.) { - mDead.push_back(entry); - } else if (entry.fraction < 0.999 || entry.ownFraction < 0.999) { - mPartial.push_back(entry); - } +/// Does the TGeo navigator's current path pass through this placement, and did it +/// stop exactly there? Compared node by node rather than as a path prefix: a +/// string prefix accepts `.../X_1` for `.../X_10`, and copy numbers 1 and 10 in +/// one mother are common enough in ALICE for that to matter. +bool tgeoPassesThrough(TGeoNavigator* nav, const std::vector& chain, bool& exact) +{ + const int depth = (int)chain.size() - 1; + const int level = nav->GetLevel(); + if (level < depth) { + return false; + } + for (int d = 0; d <= depth; ++d) { + if (nav->GetMother(level - d) != chain[d]) { + return false; } } + exact = (level == depth); + return true; +} - for (int i = 0; i < node->GetNdaughters(); ++i) { - walk(node->GetDaughter(i), here, myPath); +void sampleTask(const ReachTask& task, size_t index, int samples, TGeoNavigator* nav, ReachResult& out) +{ + TGeoVolume* volume = task.node->GetVolume(); + const bool hasDaughters = volume->GetNdaughters() > 0; + TRandom3 random(seedFor(index)); + for (int i = 0; i < samples; ++i) { + double local[3], global[3]; + if (!samplePoint(volume->GetShape(), random, local)) { + break; + } + ++out.drawn; + // nominally this volume's own material: inside its shape, inside none of its + // daughters. A leaf owns every point of its shape, so skip the walk there. + const bool own = !hasDaughters || !insideAnyDaughter(volume, local); + if (own) { + ++out.ownDrawn; + } + task.matrix.LocalToMaster(local, global); + + // FindNode() resumes from wherever the navigator currently is, so without + // this the audit asks each question from inside the very placement it is + // testing and that placement wins every genuinely ambiguous point. Two + // mutually overlapping volumes then both report themselves fully reached. + // Starting from the top makes the answer the navigator's own, and the same + // one a track crossing the region would get. + nav->CdTop(); + if (nav->FindNode(global[0], global[1], global[2]) == nullptr) { + continue; + } + bool exact = false; + if (!tgeoPassesThrough(nav, task.chain, exact)) { + continue; + } + ++out.reached; + if (own && exact) { + ++out.ownReached; // it stopped here, so the material really is this one's + } } + out.sampled = out.drawn > 0; } /// Prints the audit and returns how many placements the navigator cannot reach at all. -long reportReachability(int samples, Report& report) +long reportReachability(int samples, int jobs, Report& report) { if (samples <= 0) { return 0; } progress("reachability: asking the navigator to find every placement from inside its own shape"); - ReachAudit audit(samples); - audit.walk(gGeoManager->GetTopNode()); + ReachCollector collector; + collector.walk(gGeoManager->GetTopNode()); + auto& tasks = collector.tasks(); + + int threads = jobs > 0 ? jobs : (int)std::thread::hardware_concurrency(); + threads = std::max(1, std::min(threads, (int)tasks.size())); + // Every placement is sampled independently, so the only shared state is the + // geometry itself. ROOT serves that per thread: SetMaxThreads allocates the + // per-thread shape data (composite shapes and voxel finders cache into it) and + // each worker claims its own navigator, without which they would all drive one. + if (threads > 1) { + gGeoManager->SetMaxThreads(threads); + } + progress(form("reachability: %zu placements, %d point%s each, %d thread%s", tasks.size(), samples, + samples == 1 ? "" : "s", threads, threads == 1 ? "" : "s")); + + std::vector results(tasks.size()); + std::atomic next{0}; + auto worker = [&]() { + TGeoNavigator* nav = threads > 1 ? gGeoManager->AddNavigator() : gGeoManager->GetCurrentNavigator(); + // handed out one at a time: a placement's cost spans orders of magnitude, so a + // static split would leave most threads waiting on the few expensive ones + for (size_t i = next++; i < tasks.size(); i = next++) { + sampleTask(tasks[i], i, samples, nav, results[i]); + } + }; + if (threads > 1) { + std::vector pool; + pool.reserve(threads); + for (int i = 0; i < threads; ++i) { + pool.emplace_back(worker); + } + for (auto& thread : pool) { + thread.join(); + } + } else { + worker(); + } + + long sampled = 0, unsampleable = 0; + std::vector dead, partial; + for (size_t i = 0; i < tasks.size(); ++i) { + const auto& result = results[i]; + if (!result.sampled) { + ++unsampleable; // a sliver too thin for the rejection budget; says nothing + continue; + } + ++sampled; + TGeoNode* node = tasks[i].node; + auto* medium = node->GetVolume()->GetMedium(); + Reach entry; + entry.medium = medium != nullptr ? medium->GetName() : "(none)"; + entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; + entry.worstPath = tasks[i].path; + entry.sampled = result.drawn; + entry.ownSampled = result.ownDrawn; + entry.fraction = double(result.reached) / result.drawn; + entry.ownFraction = result.ownDrawn > 0 ? double(result.ownReached) / result.ownDrawn : 1.; + // Either number can fail on its own. A mother almost entirely filled by its + // daughters keeps a high reached fraction while the sliver of its own medium + // is taken by a foreign volume, and that sliver is the material that + // disappears -- so classify on whichever of the two is worse. + if (entry.fraction == 0.) { + dead.push_back(entry); + } else if (entry.fraction < 0.999 || entry.ownFraction < 0.999) { + partial.push_back(entry); + } + } report(form("reachability: %ld node objects visited, %ld sampled, %ld too thin to sample", - audit.nodesVisited(), audit.nodesSampled(), audit.nodesUnsampleable())); - report(form(" %ld placements the navigator never reaches, %zu it reaches only in part", - (long)audit.dead().size(), audit.partial().size())); + collector.nodesVisited(), sampled, unsampleable)); + report(form(" %ld placements the navigator never reaches, %zu it reaches only in part", (long)dead.size(), + partial.size())); // worst first: with hundreds of small overlaps the walk order is not a ranking - auto partial = audit.partial(); std::sort(partial.begin(), partial.end(), [](const Reach& a, const Reach& b) { return std::min(a.fraction, a.ownFraction) < std::min(b.fraction, b.ownFraction); }); - if (!audit.dead().empty()) { + if (!dead.empty()) { report(" unreachable -- these carry no material and produce no hits:"); report(form(" %-12s %-18s %10s %s", "mother", "medium", "sampled", "path")); - for (const auto& entry : audit.dead()) { + for (const auto& entry : dead) { report(form(" %-12s %-18s %10ld %s", entry.mother.c_str(), entry.medium.c_str(), entry.sampled, entry.worstPath.c_str())); } @@ -918,8 +1009,9 @@ long reportReachability(int samples, Report& report) report(form(" ... and %zu more, all above %.1f%%", partial.size() - 20, 100. * std::min(partial[19].fraction, partial[19].ownFraction))); } + report(""); - return (long)audit.dead().size(); + return (long)dead.size(); } // --------------------------------------------------------------------------- @@ -1739,6 +1831,7 @@ struct Options { std::vector thresholdsGauss; double margin = 5.0; int reachSamples = 32; + int reachJobs = 0; bool reachabilityOnly = false; std::string outputPrefix = "geometry-doctor"; }; @@ -1772,6 +1865,8 @@ int main(int argc, char** argv) ("reachability-samples", bpo::value(&options.reachSamples)->default_value(1000), // "points drawn inside each placement for the reachability audit; 0 disables it. Below a few " // "hundred the audit reports genuine placements as partially shadowed") // + ("reachability-jobs", bpo::value(&options.reachJobs)->default_value(0), // + "threads for the reachability audit; 0 uses every core. The answer does not depend on it") // ("reachability-only", bpo::bool_switch(&options.reachabilityOnly), // "run only the reachability audit, which needs no magnetic field"); @@ -1807,7 +1902,7 @@ int main(int argc, char** argv) report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), gGeoManager->GetListOfMedia()->GetEntries())); report(""); - const long dead = reportReachability(options.reachSamples, report); + const long dead = reportReachability(options.reachSamples, options.reachJobs, report); const std::string reportPath = options.outputPrefix + "-report.txt"; report("wrote " + reportPath); report.write(reportPath); @@ -1924,7 +2019,7 @@ int main(int argc, char** argv) report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), gGeoManager->GetListOfMedia()->GetEntries())); report(""); - reportReachability(options.reachSamples, report); + reportReachability(options.reachSamples, options.reachJobs, report); Doctor doctor(field, support); doctor.walk(gGeoManager->GetTopNode()); From 32e6a701cdd339fb65ab19558e133dc4e9e0fa99 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 4 Sep 2026 09:23:05 +0200 Subject: [PATCH 3/3] Cross-check the geometry doctor's reachability audit against VecGeom This lets the reachability audit ask a second navigator the same question, so that a placement's shadowing can be told apart from one navigator's arbitrary choice inside an overlap. - GeometryManager gains ensureVecGeomWorld() and vecGeomLocate(), which return the TGeo nodes of the path VecGeom locates a point in. Neither exposes a VecGeom type, and both answer at runtime, so a caller that does not see the private O2_WITH_VECGEOM define can still use them. - The doctor takes --navigator tgeo|vecgeom|both. 'both' reports the placements where the two engines disagree about who owns a point. - The converter flattens assemblies, so a VecGeom path is the TGeo path with its assembly levels removed and a flattened node is renamed _assemblyinternalcount_. The comparison accounts for both. - On the full ALICE geometry at 1000 samples, 575 of 27647 placements disagree. Both navigators independently find the same five unreachable MCH and beam pipe placements; they differ on the HMPID absorbers, which TGeo gives 9.7 % of their own material and VecGeom gives all of it. Co-Authored-By: Claude Opus 5 --- .../include/DetectorsBase/GeometryManager.h | 13 ++ Detectors/Base/src/GeometryManager.cxx | 45 ++++ run/CMakeLists.txt | 3 +- run/o2sim_geometry_doctor.cxx | 207 +++++++++++++++--- 4 files changed, 236 insertions(+), 32 deletions(-) diff --git a/Detectors/Base/include/DetectorsBase/GeometryManager.h b/Detectors/Base/include/DetectorsBase/GeometryManager.h index f105d137c8742..93f3931e203d4 100644 --- a/Detectors/Base/include/DetectorsBase/GeometryManager.h +++ b/Detectors/Base/include/DetectorsBase/GeometryManager.h @@ -27,6 +27,7 @@ #include "MathUtils/Cartesian.h" #include "DetectorsBase/MatCell.h" #include +#include class TGeoHMatrix; // lines 11-11 class TGeoManager; // lines 9-9 class TGeoNavigator; @@ -138,6 +139,18 @@ class GeometryManager : public TObject static constexpr bool isVecGeomAvailable() { return false; } #endif + /// Builds the VecGeom world from the currently loaded TGeo geometry, once per process, + /// and reports whether a VecGeom navigator is available at all. Unlike + /// isVecGeomAvailable() this is a runtime answer, so a caller outside this library -- + /// which does not see the private O2_WITH_VECGEOM define -- can still ask. + static bool ensureVecGeomWorld(); + + /// The VecGeom navigator's answer for a point: fills \p chain with the TGeo nodes of the + /// located path, top node first. False when this build has no VecGeom backend or the + /// point lies outside the world. Assemblies are flattened in the VecGeom geometry, so + /// the chain is shorter than the TGeo path through the same point. + static bool vecGeomLocate(double x, double y, double z, std::vector& chain); + private: /// Default constructor GeometryManager() = default; diff --git a/Detectors/Base/src/GeometryManager.cxx b/Detectors/Base/src/GeometryManager.cxx index 225f21d8239a1..69c5a8d25498b 100644 --- a/Detectors/Base/src/GeometryManager.cxx +++ b/Detectors/Base/src/GeometryManager.cxx @@ -678,3 +678,48 @@ o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, f } #endif // O2_WITH_VECGEOM + +//_____________________________________________________________________________________ +bool GeometryManager::ensureVecGeomWorld() +{ +#ifdef O2_WITH_VECGEOM + ensureVecGeomWorldBuilt(); + return true; +#else + return false; +#endif +} + +//_____________________________________________________________________________________ +bool GeometryManager::vecGeomLocate(double x, double y, double z, std::vector& chain) +{ + chain.clear(); +#ifdef O2_WITH_VECGEOM + ensureVecGeomWorldBuilt(); + // One state per thread, as for the material budget above. + thread_local vecgeom::NavigationState* state = + vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth()); + state->Clear(); + const vecgeom::Vector3D point(x, y, z); + if (vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), point, *state, true) == + nullptr) { + return false; + } + auto const& converter = tgeo2vecgeom::RootGeoManager::Instance(); + for (int level = 0; level < (int)state->GetCurrentLevel(); ++level) { + auto const* placed = state->At(level); + auto const* node = placed != nullptr ? converter.tgeonode(placed) : nullptr; + if (node == nullptr) { + chain.clear(); + return false; + } + chain.push_back(const_cast(node)); + } + return !chain.empty(); +#else + (void)x; + (void)y; + (void)z; + return false; +#endif +} diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index db00f0465e645..0b88a6e6f68d5 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -115,7 +115,8 @@ o2_add_executable(g4-determine-unknown-pdg-properties o2_add_executable(geometry-doctor COMPONENT_NAME sim SOURCES o2sim_geometry_doctor.cxx - PUBLIC_LINK_LIBRARIES O2::Field ROOT::Geom Boost::program_options + PUBLIC_LINK_LIBRARIES O2::Field O2::DetectorsBase ROOT::Geom + Boost::program_options nlohmann_json::nlohmann_json) o2_add_executable(mctracks-proxy diff --git a/run/o2sim_geometry_doctor.cxx b/run/o2sim_geometry_doctor.cxx index 87d41a4f83167..980d7c711411d 100644 --- a/run/o2sim_geometry_doctor.cxx +++ b/run/o2sim_geometry_doctor.cxx @@ -41,6 +41,7 @@ /// sensitivity are GSTMED parameters 1 and 0 of the TGeoMedium, and are recovered /// from the file itself. +#include "DetectorsBase/GeometryManager.h" #include "Field/MagneticField.h" #include @@ -737,10 +738,21 @@ struct Reach { double fraction = 1.; long ownSampled = 0; ///< points that are nominally this volume's own material double ownFraction = 1.; ///< of those, the share the navigator actually gives it + double vgFraction = 1.; ///< the same two numbers from the VecGeom navigator, if asked + double vgOwnFraction = 1.; + long disagreed = 0; ///< points the two navigators answer differently }; constexpr int kReachRejectionTries = 400; +/// Which navigator answers "what is at this point". They are independent +/// implementations of the same question, and TGeo's answer inside an overlap is +/// arbitrary -- so where the two disagree, the overlap is not only real but its +/// resolution depends on which engine the simulation ran with. +enum class Navigator { TGeo, + VecGeom, + Both }; + /// Defined with the placement table below. Deliberately the same predicate the /// field classification already uses for "own material", so the two parts of this /// tool cannot disagree about what a volume's own material is. @@ -754,12 +766,14 @@ struct ReachTask { TGeoNode* node = nullptr; TGeoHMatrix matrix; std::string path; - std::vector chain; ///< top node first, this node last + std::vector chain; ///< top node first, this node last + std::vector flatChain; ///< the same, with the assembly levels removed }; struct ReachResult { bool sampled = false; long drawn = 0, reached = 0, ownDrawn = 0, ownReached = 0; + long vgReached = 0, vgOwnReached = 0, disagreed = 0; }; /// Collects one representative placement per node object. This half of the audit @@ -768,19 +782,20 @@ struct ReachResult { class ReachCollector { public: - void walk(TGeoNode* node) { walk(node, TGeoHMatrix(), "", {}); } + void walk(TGeoNode* node) { walk(node, TGeoHMatrix(), "", {}, {}); } std::vector& tasks() { return mTasks; } long nodesVisited() const { return mVisited; } private: - void walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path, std::vector chain); + void walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path, std::vector chain, + std::vector flatChain); std::set mSeen; std::vector mTasks; long mVisited = 0; }; void ReachCollector::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path, - std::vector chain) + std::vector chain, std::vector flatChain) { if (!mSeen.insert(node).second) { return; // this node object, and therefore its whole subtree, is already covered @@ -793,15 +808,17 @@ void ReachCollector::walk(TGeoNode* node, const TGeoHMatrix& parent, const std:: // an assembly is expanded away at closure, so FindNode never returns one if (!node->GetVolume()->IsAssembly()) { + flatChain.push_back(node); ReachTask task; task.node = node; task.matrix = here; task.path = myPath; task.chain = chain; + task.flatChain = flatChain; mTasks.push_back(std::move(task)); } for (int i = 0; i < node->GetNdaughters(); ++i) { - walk(node->GetDaughter(i), here, myPath, chain); + walk(node->GetDaughter(i), here, myPath, chain, flatChain); } } @@ -861,10 +878,52 @@ bool tgeoPassesThrough(TGeoNavigator* nav, const std::vector& chain, return true; } -void sampleTask(const ReachTask& task, size_t index, int samples, TGeoNavigator* nav, ReachResult& out) +/// The converter flattens assemblies the way Geant4 does, and a flattened daughter +/// becomes a *new* TGeoNode -- same volume, same place, name +/// `_assemblyinternalcount_`. So a located node is this one either by +/// pointer, or by carrying its volume and its name under that suffix. The rule +/// cannot separate two placements of one volume that share a node name, which the +/// MFT half cone does; there the VecGeom answer is the weaker of the two. +bool sameFlattenedNode(TGeoNode* located, TGeoNode* wanted) +{ + if (located == wanted) { + return true; + } + if (located->GetVolume() != wanted->GetVolume()) { + return false; + } + static const std::string kFlattened = "_assemblyinternalcount_"; + const std::string name = located->GetName(), want = wanted->GetName(); + return name.size() > want.size() + kFlattened.size() && name.compare(0, want.size(), want) == 0 && + name.compare(want.size(), kFlattened.size(), kFlattened) == 0; +} + +/// The same question of the VecGeom path. Assemblies carry no material and are gone +/// from the VecGeom geometry, so the comparison is against the TGeo path with its +/// assembly levels removed -- which is what the two navigators genuinely have in +/// common. +bool vecGeomPassesThrough(const std::vector& located, const std::vector& flatChain, bool& exact) +{ + const int depth = (int)flatChain.size() - 1; + if ((int)located.size() - 1 < depth) { + return false; + } + for (int d = 0; d <= depth; ++d) { + if (!sameFlattenedNode(located[d], flatChain[d])) { + return false; + } + } + exact = ((int)located.size() - 1 == depth); + return true; +} + +void sampleTask(const ReachTask& task, size_t index, int samples, Navigator backend, TGeoNavigator* nav, + std::vector& located, ReachResult& out) { TGeoVolume* volume = task.node->GetVolume(); const bool hasDaughters = volume->GetNdaughters() > 0; + const bool wantTGeo = backend != Navigator::VecGeom; + const bool wantVecGeom = backend != Navigator::TGeo; TRandom3 random(seedFor(index)); for (int i = 0; i < samples; ++i) { double local[3], global[3]; @@ -880,30 +939,46 @@ void sampleTask(const ReachTask& task, size_t index, int samples, TGeoNavigator* } task.matrix.LocalToMaster(local, global); - // FindNode() resumes from wherever the navigator currently is, so without - // this the audit asks each question from inside the very placement it is - // testing and that placement wins every genuinely ambiguous point. Two - // mutually overlapping volumes then both report themselves fully reached. - // Starting from the top makes the answer the navigator's own, and the same - // one a track crossing the region would get. - nav->CdTop(); - if (nav->FindNode(global[0], global[1], global[2]) == nullptr) { - continue; - } - bool exact = false; - if (!tgeoPassesThrough(nav, task.chain, exact)) { - continue; + bool tgeoThrough = false, tgeoExact = false; + if (wantTGeo) { + // FindNode() resumes from wherever the navigator currently is, so without + // this the audit asks each question from inside the very placement it is + // testing and that placement wins every genuinely ambiguous point. Two + // mutually overlapping volumes then both report themselves fully reached. + // Starting from the top makes the answer the navigator's own, and the same + // one a track crossing the region would get. + nav->CdTop(); + if (nav->FindNode(global[0], global[1], global[2]) != nullptr) { + tgeoThrough = tgeoPassesThrough(nav, task.chain, tgeoExact); + } + if (tgeoThrough) { + ++out.reached; + if (own && tgeoExact) { + ++out.ownReached; // it stopped here, so the material really is this one's + } + } } - ++out.reached; - if (own && exact) { - ++out.ownReached; // it stopped here, so the material really is this one's + if (wantVecGeom) { + bool vgThrough = false, vgExact = false; + if (o2::base::GeometryManager::vecGeomLocate(global[0], global[1], global[2], located)) { + vgThrough = vecGeomPassesThrough(located, task.flatChain, vgExact); + } + if (vgThrough) { + ++out.vgReached; + if (own && vgExact) { + ++out.vgOwnReached; + } + } + if (backend == Navigator::Both && vgThrough != tgeoThrough) { + ++out.disagreed; + } } } out.sampled = out.drawn > 0; } /// Prints the audit and returns how many placements the navigator cannot reach at all. -long reportReachability(int samples, int jobs, Report& report) +long reportReachability(int samples, int jobs, Navigator backend, Report& report) { if (samples <= 0) { return 0; @@ -913,6 +988,11 @@ long reportReachability(int samples, int jobs, Report& report) collector.walk(gGeoManager->GetTopNode()); auto& tasks = collector.tasks(); + if (backend != Navigator::TGeo && !o2::base::GeometryManager::ensureVecGeomWorld()) { + report(" VecGeom backend requested but this build of O2 has none; falling back to TGeo"); + backend = Navigator::TGeo; + } + int threads = jobs > 0 ? jobs : (int)std::thread::hardware_concurrency(); threads = std::max(1, std::min(threads, (int)tasks.size())); // Every placement is sampled independently, so the only shared state is the @@ -929,10 +1009,11 @@ long reportReachability(int samples, int jobs, Report& report) std::atomic next{0}; auto worker = [&]() { TGeoNavigator* nav = threads > 1 ? gGeoManager->AddNavigator() : gGeoManager->GetCurrentNavigator(); + std::vector located; // handed out one at a time: a placement's cost spans orders of magnitude, so a // static split would leave most threads waiting on the few expensive ones for (size_t i = next++; i < tasks.size(); i = next++) { - sampleTask(tasks[i], i, samples, nav, results[i]); + sampleTask(tasks[i], i, samples, backend, nav, located, results[i]); } }; if (threads > 1) { @@ -948,7 +1029,7 @@ long reportReachability(int samples, int jobs, Report& report) worker(); } - long sampled = 0, unsampleable = 0; + long sampled = 0, unsampleable = 0, disagreeing = 0; std::vector dead, partial; for (size_t i = 0; i < tasks.size(); ++i) { const auto& result = results[i]; @@ -965,8 +1046,17 @@ long reportReachability(int samples, int jobs, Report& report) entry.worstPath = tasks[i].path; entry.sampled = result.drawn; entry.ownSampled = result.ownDrawn; - entry.fraction = double(result.reached) / result.drawn; - entry.ownFraction = result.ownDrawn > 0 ? double(result.ownReached) / result.ownDrawn : 1.; + entry.disagreed = result.disagreed; + const bool primaryIsVecGeom = backend == Navigator::VecGeom; + entry.fraction = double(primaryIsVecGeom ? result.vgReached : result.reached) / result.drawn; + entry.ownFraction = result.ownDrawn > 0 + ? double(primaryIsVecGeom ? result.vgOwnReached : result.ownReached) / result.ownDrawn + : 1.; + entry.vgFraction = double(result.vgReached) / result.drawn; + entry.vgOwnFraction = result.ownDrawn > 0 ? double(result.vgOwnReached) / result.ownDrawn : 1.; + if (result.disagreed > 0) { + ++disagreeing; + } // Either number can fail on its own. A mother almost entirely filled by its // daughters keeps a high reached fraction while the sliver of its own medium // is taken by a foreign volume, and that sliver is the material that @@ -978,8 +1068,9 @@ long reportReachability(int samples, int jobs, Report& report) } } - report(form("reachability: %ld node objects visited, %ld sampled, %ld too thin to sample", - collector.nodesVisited(), sampled, unsampleable)); + const char* named = backend == Navigator::TGeo ? "TGeo" : (backend == Navigator::VecGeom ? "VecGeom" : "TGeo, cross-checked against VecGeom"); + report(form("reachability: %ld node objects visited, %ld sampled, %ld too thin to sample (navigator: %s)", + collector.nodesVisited(), sampled, unsampleable, named)); report(form(" %ld placements the navigator never reaches, %zu it reaches only in part", (long)dead.size(), partial.size())); // worst first: with hundreds of small overlaps the walk order is not a ranking @@ -1010,6 +1101,46 @@ long reportReachability(int samples, int jobs, Report& report) 100. * std::min(partial[19].fraction, partial[19].ownFraction))); } + if (backend == Navigator::Both) { + report(""); + if (disagreeing == 0) { + report(" TGeo and VecGeom agree on every point sampled."); + } else { + report(form(" %ld placements where the two navigators disagree about who owns a point.", disagreeing)); + report(" A disagreement is a real overlap whose resolution depends on the engine, so the"); + report(" material a track sees there is not a property of the geometry alone:"); + report(form(" %-12s %-18s %8s %8s %8s %s", "mother", "medium", "differ", "TGeo", "VecGeom", "path")); + std::vector conflicts; + for (size_t i = 0; i < tasks.size(); ++i) { + if (results[i].disagreed == 0 || !results[i].sampled) { + continue; + } + Reach entry; + TGeoNode* node = tasks[i].node; + auto* medium = node->GetVolume()->GetMedium(); + entry.medium = medium != nullptr ? medium->GetName() : "(none)"; + entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; + entry.worstPath = tasks[i].path; + entry.sampled = results[i].drawn; + entry.disagreed = results[i].disagreed; + entry.fraction = double(results[i].reached) / results[i].drawn; + entry.vgFraction = double(results[i].vgReached) / results[i].drawn; + conflicts.push_back(entry); + } + std::sort(conflicts.begin(), conflicts.end(), [](const Reach& a, const Reach& b) { + return double(a.disagreed) / a.sampled > double(b.disagreed) / b.sampled; + }); + for (size_t i = 0; i < conflicts.size() && i < 20; ++i) { + const auto& entry = conflicts[i]; + report(form(" %-12s %-18s %7.1f%% %7.1f%% %7.1f%% %s", entry.mother.c_str(), entry.medium.c_str(), + 100. * entry.disagreed / entry.sampled, 100. * entry.fraction, 100. * entry.vgFraction, + entry.worstPath.c_str())); + } + if (conflicts.size() > 20) { + report(form(" ... and %zu more", conflicts.size() - 20)); + } + } + } report(""); return (long)dead.size(); } @@ -1832,6 +1963,7 @@ struct Options { double margin = 5.0; int reachSamples = 32; int reachJobs = 0; + std::string navigator = "tgeo"; bool reachabilityOnly = false; std::string outputPrefix = "geometry-doctor"; }; @@ -1867,6 +1999,9 @@ int main(int argc, char** argv) "hundred the audit reports genuine placements as partially shadowed") // ("reachability-jobs", bpo::value(&options.reachJobs)->default_value(0), // "threads for the reachability audit; 0 uses every core. The answer does not depend on it") // + ("navigator", bpo::value(&options.navigator)->default_value("tgeo"), // + "which navigator answers 'what is at this point': tgeo, vecgeom, or both. 'both' reports where " // + "they disagree, which is where a real overlap is resolved differently by the two engines") // ("reachability-only", bpo::bool_switch(&options.reachabilityOnly), // "run only the reachability audit, which needs no magnetic field"); @@ -1884,6 +2019,16 @@ int main(int argc, char** argv) return 1; } + Navigator navigator = Navigator::TGeo; + if (options.navigator == "vecgeom") { + navigator = Navigator::VecGeom; + } else if (options.navigator == "both") { + navigator = Navigator::Both; + } else if (options.navigator != "tgeo") { + std::cerr << "error: --navigator takes tgeo, vecgeom or both\n"; + return 1; + } + const bool haveFieldFile = arguments.count("field-file") != 0u; const bool haveFieldCurrent = arguments.count("field-current") != 0u; @@ -1902,7 +2047,7 @@ int main(int argc, char** argv) report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), gGeoManager->GetListOfMedia()->GetEntries())); report(""); - const long dead = reportReachability(options.reachSamples, options.reachJobs, report); + const long dead = reportReachability(options.reachSamples, options.reachJobs, navigator, report); const std::string reportPath = options.outputPrefix + "-report.txt"; report("wrote " + reportPath); report.write(reportPath); @@ -2019,7 +2164,7 @@ int main(int argc, char** argv) report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), gGeoManager->GetListOfMedia()->GetEntries())); report(""); - reportReachability(options.reachSamples, options.reachJobs, report); + reportReachability(options.reachSamples, options.reachJobs, navigator, report); Doctor doctor(field, support); doctor.walk(gGeoManager->GetTopNode());