@@ -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.
721731struct 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
727739constexpr 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+
729756class 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