Skip to content

Commit d4c9de5

Browse files
sawenzelclaude
andcommitted
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 <noreply@anthropic.com>
1 parent 9d5742d commit d4c9de5

1 file changed

Lines changed: 84 additions & 24 deletions

File tree

run/o2sim_geometry_doctor.cxx

Lines changed: 84 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -718,14 +718,41 @@ bool supportFromJson(const json& in, Support& support)
718718
/// common defect lives -- a daughter outside its mother is a property of the node,
719719
/// not of the path that reaches it. A node whose mother is itself placed many
720720
/// times is therefore sampled once, in the first of those placements.
721+
///
722+
/// Two questions are asked of every point, and they are not the same question.
723+
/// *Reachability* asks whether the navigator's path passes through this placement
724+
/// at all, so a point that lands in one of its own daughters counts. *Self
725+
/// material* asks the stronger question: for a point that is nominally this
726+
/// volume's own material -- inside its shape and inside none of its daughters --
727+
/// FindNode() must return exactly this path, not a prefix of it and not something
728+
/// else. A mother whose own medium is entirely taken by an overlapping foreign
729+
/// volume is still "reached" through its daughters, and only the second question
730+
/// sees that its material is gone.
721731
struct Reach {
722732
std::string medium, mother, worstPath;
723733
long sampled = 0;
724734
double fraction = 1.;
735+
long ownSampled = 0; ///< points that are nominally this volume's own material
736+
double ownFraction = 1.; ///< of those, the share the navigator actually gives it
725737
};
726738

727739
constexpr int kReachRejectionTries = 400;
728740

741+
/// `found` passes through `path` -- a prefix match that must end on a path
742+
/// separator. Without the boundary check `.../X_1` matches `.../X_10`, and copy
743+
/// numbers 1 and 10 in one mother are common enough in ALICE that the check would
744+
/// silently accept a point the navigator gave to a different sibling.
745+
inline bool passesThrough(const std::string& found, const std::string& path)
746+
{
747+
return found.compare(0, path.size(), path) == 0 &&
748+
(found.size() == path.size() || found[path.size()] == '/');
749+
}
750+
751+
/// Defined with the placement table below. Deliberately the same predicate the
752+
/// field classification already uses for "own material", so the two parts of this
753+
/// tool cannot disagree about what a volume's own material is.
754+
bool insideAnyDaughter(TGeoVolume* volume, const double* local);
755+
729756
class ReachAudit
730757
{
731758
public:
@@ -785,37 +812,61 @@ void ReachAudit::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::stri
785812

786813
// an assembly is expanded away at closure, so FindNode never returns one
787814
if (!volume->IsAssembly()) {
788-
int drawn = 0, reached = 0;
815+
int drawn = 0, reached = 0, ownDrawn = 0, ownReached = 0;
816+
const bool hasDaughters = volume->GetNdaughters() > 0;
789817
for (int i = 0; i < mSamples; ++i) {
790818
double local[3], global[3];
791819
if (!samplePoint(volume->GetShape(), local)) {
792820
break;
793821
}
794822
++drawn;
823+
// nominally this volume's own material: inside its shape, inside none of its
824+
// daughters. A leaf owns every point of its shape, so skip the walk there.
825+
const bool own = !hasDaughters || !insideAnyDaughter(volume, local);
826+
if (own) {
827+
++ownDrawn;
828+
}
795829
here.LocalToMaster(local, global);
830+
// FindNode() resumes from wherever the navigator currently is, so without
831+
// this the audit asks each question from inside the very placement it is
832+
// testing and that placement wins every genuinely ambiguous point. Two
833+
// mutually overlapping volumes then both report themselves fully reached.
834+
// Starting from the top makes the answer the navigator's own, and the same
835+
// one a track crossing the region would get.
836+
gGeoManager->CdTop();
796837
if (gGeoManager->FindNode(global[0], global[1], global[2]) == nullptr) {
797838
continue;
798839
}
799840
const std::string found = gGeoManager->GetPath();
800-
// reached if the navigator's own path passes through this placement
801-
if (found.rfind(myPath, 0) == 0) {
802-
++reached;
841+
if (!passesThrough(found, myPath)) {
842+
continue;
843+
}
844+
++reached; // the navigator's own path passes through this placement
845+
if (own && found.size() == myPath.size()) {
846+
++ownReached; // ...and it stopped here, so the material really is this one's
803847
}
804848
}
805849
if (drawn == 0) {
806850
++mUnsampleable; // a sliver too thin for the rejection budget; says nothing
807851
} else {
808852
++mSampled;
809-
const double fraction = double(reached) / drawn;
810-
if (fraction < 0.999) {
811-
auto* medium = volume->GetMedium();
812-
Reach entry;
813-
entry.medium = medium != nullptr ? medium->GetName() : "(none)";
814-
entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-";
815-
entry.worstPath = myPath;
816-
entry.sampled = drawn;
817-
entry.fraction = fraction;
818-
(fraction == 0. ? mDead : mPartial).push_back(entry);
853+
auto* medium = volume->GetMedium();
854+
Reach entry;
855+
entry.medium = medium != nullptr ? medium->GetName() : "(none)";
856+
entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-";
857+
entry.worstPath = myPath;
858+
entry.sampled = drawn;
859+
entry.fraction = double(reached) / drawn;
860+
entry.ownSampled = ownDrawn;
861+
entry.ownFraction = ownDrawn > 0 ? double(ownReached) / ownDrawn : 1.;
862+
// Either number can fail on its own. A mother almost entirely filled by its
863+
// daughters keeps a high reached fraction while the sliver of its own medium
864+
// is taken by a foreign volume, and that sliver is the material that
865+
// disappears -- so classify on whichever of the two is worse.
866+
if (entry.fraction == 0.) {
867+
mDead.push_back(entry);
868+
} else if (entry.fraction < 0.999 || entry.ownFraction < 0.999) {
869+
mPartial.push_back(entry);
819870
}
820871
}
821872
}
@@ -839,6 +890,11 @@ long reportReachability(int samples, Report& report)
839890
audit.nodesVisited(), audit.nodesSampled(), audit.nodesUnsampleable()));
840891
report(form(" %ld placements the navigator never reaches, %zu it reaches only in part",
841892
(long)audit.dead().size(), audit.partial().size()));
893+
// worst first: with hundreds of small overlaps the walk order is not a ranking
894+
auto partial = audit.partial();
895+
std::sort(partial.begin(), partial.end(), [](const Reach& a, const Reach& b) {
896+
return std::min(a.fraction, a.ownFraction) < std::min(b.fraction, b.ownFraction);
897+
});
842898
if (!audit.dead().empty()) {
843899
report(" unreachable -- these carry no material and produce no hits:");
844900
report(form(" %-12s %-18s %10s %s", "mother", "medium", "sampled", "path"));
@@ -847,17 +903,20 @@ long reportReachability(int samples, Report& report)
847903
entry.worstPath.c_str()));
848904
}
849905
}
850-
for (size_t i = 0; i < audit.partial().size() && i < 20; ++i) {
851-
const auto& entry = audit.partial()[i];
906+
for (size_t i = 0; i < partial.size() && i < 20; ++i) {
907+
const auto& entry = partial[i];
852908
if (i == 0) {
853-
report(" partially shadowed -- an overlapping sibling or an extruding placement:");
854-
report(form(" %-12s %-18s %8s %s", "mother", "medium", "reached", "path"));
909+
report(" partially shadowed -- an overlapping sibling or an extruding placement.");
910+
report(" 'reached' is how much of the placement the navigator enters at all; 'own kept'");
911+
report(" how much of the medium this volume was given to carry survives as its own:");
912+
report(form(" %-12s %-18s %8s %9s %s", "mother", "medium", "reached", "own kept", "path"));
855913
}
856-
report(form(" %-12s %-18s %7.1f%% %s", entry.mother.c_str(), entry.medium.c_str(),
857-
100. * entry.fraction, entry.worstPath.c_str()));
914+
report(form(" %-12s %-18s %7.1f%% %8.1f%% %s", entry.mother.c_str(), entry.medium.c_str(),
915+
100. * entry.fraction, 100. * entry.ownFraction, entry.worstPath.c_str()));
858916
}
859-
if (audit.partial().size() > 20) {
860-
report(form(" ... and %zu more", audit.partial().size() - 20));
917+
if (partial.size() > 20) {
918+
report(form(" ... and %zu more, all above %.1f%%", partial.size() - 20,
919+
100. * std::min(partial[19].fraction, partial[19].ownFraction)));
861920
}
862921
report("");
863922
return (long)audit.dead().size();
@@ -1710,8 +1769,9 @@ int main(int argc, char** argv)
17101769
"prefix for the report, the proposals and the placement table") //
17111770
("verify-anchors", bpo::value<std::string>(&options.anchorFile), //
17121771
"check the classification against known-good volumes listed in this JSON file") //
1713-
("reachability-samples", bpo::value<int>(&options.reachSamples)->default_value(32), //
1714-
"points drawn inside each placement for the reachability audit; 0 disables it") //
1772+
("reachability-samples", bpo::value<int>(&options.reachSamples)->default_value(1000), //
1773+
"points drawn inside each placement for the reachability audit; 0 disables it. Below a few " //
1774+
"hundred the audit reports genuine placements as partially shadowed") //
17151775
("reachability-only", bpo::bool_switch(&options.reachabilityOnly), //
17161776
"run only the reachability audit, which needs no magnetic field");
17171777

0 commit comments

Comments
 (0)