diff --git a/geom/geom/src/TGeoNavigator.cxx b/geom/geom/src/TGeoNavigator.cxx index 0c6a9781110de..ceca7424ee7c1 100644 --- a/geom/geom/src/TGeoNavigator.cxx +++ b/geom/geom/src/TGeoNavigator.cxx @@ -930,6 +930,39 @@ TGeoNode *TGeoNavigator::FindNextBoundary(Double_t stepmax, const char *path, Bo } matrix->MasterToLocal(fPoint, dpt); matrix->MasterToLocalVect(fDirection, dvec); + // If the node below this mother is a MANY node, its overlap + // candidates have to be checked even when the current node is + // one of its ONLY descendants. + if (ovlp && !offset) { + Int_t *ovlps = currentnode->GetOverlaps(novlps); + for (Int_t i = 0; i < novlps; i++) { + current = mothernode->GetVolume()->GetNode(ovlps[i]); + if (current->IsOverlapping()) + continue; + current->cd(); + current->MasterToLocal(dpt, mothpt); + current->MasterToLocalVect(dvec, vecpt); + snext = current->GetVolume()->GetShape()->DistFromOutside(mothpt, vecpt, iact, fStep, &safe); + if (snext < fStep - gTolerance) { + if (computeGlobal) { + fCurrentMatrix->CopyFrom(matrix); + fCurrentMatrix->Multiply(current->GetMatrix()); + } + fIsStepExiting = kTRUE; + fIsStepEntering = kFALSE; + fStep = snext; + fNextNode = current; + fNextDaughterIndex = -3; + PushPath(); + Int_t iup = up; + while (iup--) + CdUp(); + CdDown(ovlps[i]); + DoBackupState(); + PopPath(); + } + } + } snext = TGeoShape::Big(); if (!mothernode->GetVolume()->IsAssembly()) snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt, dvec, iact, fStep); @@ -1446,6 +1479,53 @@ TGeoNode *TGeoNavigator::FindNextBoundaryAndStep(Double_t stepmax, Bool_t compsa matrix = GetMotherMatrix(up); matrix->MasterToLocal(fPoint, dpt); matrix->MasterToLocalVect(fDirection, dvec); + /* + * The live leaf can be ONLY while currentnode is a MANY ancestor. + * Check its ONLY overlap candidates in mothernode. If one wins, + * commit the climb to mothernode and restart because the path, + * MANY count, and coordinate frame have changed. + */ + Bool_t restart = kFALSE; + if (ovlp && !offset) { + Int_t *ovlps = currentnode->GetOverlaps(novlps); + Int_t icandidate = -1; + TGeoNode *candidate = nullptr; + for (Int_t i = 0; i < novlps; i++) { + current = mothernode->GetVolume()->GetNode(ovlps[i]); + if (current->IsOverlapping()) + continue; + current->cd(); + current->MasterToLocal(dpt, mothpt); + current->MasterToLocalVect(dvec, vecpt); + snext = current->GetVolume()->GetShape()->DistFromOutside(mothpt, vecpt, iact, fStep); + if (snext < fStep - gTolerance) { + fCurrentMatrix->CopyFrom(matrix); + fCurrentMatrix->Multiply(current->GetMatrix()); + fIsStepEntering = kFALSE; + fIsStepExiting = kTRUE; + fStep = snext; + fNextNode = current; + icandidate = ovlps[i]; + candidate = current; + } + } + if (icandidate >= 0) { + icrossed = icandidate; + current = candidate; + Int_t iup = up; + while (iup--) + CdUp(); + PopDummy(); + PushPath(fLevel + 1); + nmany = fNmany; + up = 1; + currentnode = fCurrentNode; + ovlp = currentnode->IsOverlapping(); + restart = kTRUE; + } + if (restart) + continue; + } snext = TGeoShape::Big(); if (!mothernode->GetVolume()->IsAssembly()) snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt, dvec, iact, fStep); diff --git a/geom/test/CMakeLists.txt b/geom/test/CMakeLists.txt index afcf07a6fa6a6..c52ed9610d743 100644 --- a/geom/test/CMakeLists.txt +++ b/geom/test/CMakeLists.txt @@ -21,6 +21,10 @@ ROOT_ADD_GTEST(tessellated test_tessellated.cxx LIBRARIES Geom) +ROOT_ADD_GTEST(overlap_navigation + test_overlap_navigation.cxx + LIBRARIES Geom) + if(imt) ROOT_ADD_GTEST(manager_lifetime test_manager_lifetime.cxx diff --git a/geom/test/test_overlap_navigation.cxx b/geom/test/test_overlap_navigation.cxx new file mode 100644 index 0000000000000..88ba3c095a36c --- /dev/null +++ b/geom/test/test_overlap_navigation.cxx @@ -0,0 +1,68 @@ +#include + +#include +#include +#include +#include +#include + +#include + +TEST(TGeoNavigator, ManyOverlapIsResolvedFromOnlyDaughter) +{ + auto geom = std::make_unique("many_overlap", "MANY overlap navigation test"); + + auto *matAir = new TGeoMaterial("Air"); + auto *matAl = new TGeoMaterial("Al"); + auto *matPb = new TGeoMaterial("Pb"); + auto *air = new TGeoMedium("air", 1, matAir); + auto *al = new TGeoMedium("al", 2, matAl); + auto *pb = new TGeoMedium("pb", 3, matPb); + + auto *top = geom->MakeBox("TOP", air, 50., 50., 50.); + auto *mother = geom->MakeBox("M", air, 20., 20., 20.); + auto *many = geom->MakeBox("A", al, 5., 5., 5.); + auto *daughter = geom->MakeBox("D", air, 2., 2., 5.); + auto *onlySibling = geom->MakeBox("B", pb, 2., 2., 3.); + + // M contains overlapping siblings A (MANY) and B (ONLY), while A contains + // the ONLY daughter D. D covers the full region where B overlaps A: + // + // M + // |-- A (MANY) + // | `-- D (ONLY) + // `-- B (ONLY) + geom->SetTopVolume(top); + top->AddNode(mother, 1); + mother->AddNodeOverlap(many, 1); + mother->AddNode(onlySibling, 1); + many->AddNode(daughter, 1); + geom->CloseGeometry(); + + // Point location establishes the expected priority: B wins over A and D. + ASSERT_EQ(geom->FindNode(0., 0., 0.)->GetVolume(), onlySibling); + + // The first step enters D at z = -5. The non-stepping query must retain + // that current path while finding B at z = -3 as the next boundary. + geom->InitTrack(0., 0., -10., 0., 0., 1.); + ASSERT_STREQ(geom->GetCurrentNode()->GetName(), "M_1"); + + geom->FindNextBoundaryAndStep(); + EXPECT_DOUBLE_EQ(geom->GetStep(), 5.); + ASSERT_STREQ(geom->GetCurrentNode()->GetName(), "D_1"); + + auto *next = geom->FindNextBoundary(); + EXPECT_DOUBLE_EQ(geom->GetStep(), 2.); + ASSERT_NE(next, nullptr); + EXPECT_EQ(next->GetVolume(), onlySibling); + + // Repeat with the combined API: its second step must commit the sibling + // transition from D to B rather than crossing all of D. + geom->InitTrack(0., 0., -10., 0., 0., 1.); + geom->FindNextBoundaryAndStep(); + ASSERT_STREQ(geom->GetCurrentNode()->GetName(), "D_1"); + + geom->FindNextBoundaryAndStep(); + EXPECT_DOUBLE_EQ(geom->GetStep(), 2.); + EXPECT_STREQ(geom->GetCurrentNode()->GetName(), "B_1"); +} diff --git a/test/stressGeometry.cxx b/test/stressGeometry.cxx index 772098d927196..e82d5122941d1 100644 --- a/test/stressGeometry.cxx +++ b/test/stressGeometry.cxx @@ -127,39 +127,39 @@ const char *exps[NG] = {"aleph", "belle", "atlas" }; -const Int_t versions[NG] = {5, //aleph - 3, //barres - 3, //felix - 3, //phenix - 3, //chambers - 4, //p326 - 3, //bes - 3, //dubna - 3, //ganil - 3, //e907 - 4, //phobos2 - 3, //hermes - 3, //na35 - 3, //na47 - 3, //na49 - 3, //wa91 - 3, //sdc - 4, //integral - 4, //ams - 3, //brahms - 5, //gem - 4, //tesla - 3, //btev - 6, //cdf - 4, //hades2 - 4, //lhcbfull - 4, //star - 4, //sld - 4, //cms - 6, //alice3 - 4, //babar2 - 3, //belle - 6}; //atlas +const Int_t versions[NG] = {5, // aleph + 3, // barres + 3, // felix + 3, // phenix + 3, // chambers + 4, // p326 + 3, // bes + 3, // dubna + 3, // ganil + 3, // e907 + 4, // phobos2 + 3, // hermes + 3, // na35 + 3, // na47 + 3, // na49 + 3, // wa91 + 3, // sdc + 4, // integral + 4, // ams + 4, // brahms + 5, // gem + 4, // tesla + 3, // btev + 6, // cdf + 4, // hades2 + 4, // lhcbfull + 5, // star + 4, // sld + 4, // cms + 6, // alice3 + 5, // babar2 + 3, // belle + 7}; // atlas // The timings below are on my machine PIV 3GHz const Double_t cp_brun[NG] = {1.9, //aleph 0.1, //barres @@ -322,7 +322,7 @@ void ReadRef(Int_t kexp) { if (!gen_ref) fname = TString::Format("root://eospublic.cern.ch//eos/root-eos/testfiles//%s_ref_%d.root", exps[kexp],versions[kexp]); else - fname.Format("files/%s_ref_%d.root", exps[kexp],versions[kexp]); + fname = TString::Format("files/%s_ref_%d.root", exps[kexp], versions[kexp]); f = TFile::Open(fname,"CACHEREAD"); if (!f) {