u4/tests/U4SolidMakerTest2.cc

0001 /**
0002 U4SolidMakerTest2.cc
0003 ======================
0004
0005 Test voxelized G4MultiUnion navigation via the public Geant4 API.
0006
0007 Validates that ssolid::Distance (which now always uses the standard
0008 Inside/DistanceToIn/DistanceToOut path) returns correct results for
0009 G4MultiUnion solids produced by U4SolidMaker.
0010
0011 Exercises:
0012   - rays that clearly hit a constituent solid
0013   - rays that miss entirely
0014   - rays originating inside a constituent
0015   - rays along axes through multiple constituents
0016   - Inside() classification at known points
0017   - verification that Voxelize() was called (even if voxel count is 0
0018     for small numbers of solids, navigation must still work)
0019
0020 Returns 0 on success, non-zero on failure.
0021 **/
0022
0023 #include <iostream>
0024 #include <iomanip>
0025 #include <cmath>
0026 #include <cassert>
0027 #include <vector>
0028 #include <string>
0029
0030 #include "G4VSolid.hh"
0031 #include "G4MultiUnion.hh"
0032 #include "G4ThreeVector.hh"
0033 #include "geomdefs.hh"
0034
0035 #include "U4SolidMaker.hh"
0036 #include "ssolid.h"
0037 #include "OPTICKS_LOG.hh"
0038
0039 static const double kTol = 1.0e-3 ; // mm
0040
0041 struct Ray
0042 {
0043     const char* label ;
0044     G4ThreeVector pos ;
0045     G4ThreeVector dir ;
0046     double        expected_t ;   // kInfinity  → expect a miss
0047     double        tolerance ;
0048 };
0049
0050 struct InsidePoint
0051 {
0052     const char* label ;
0053     G4ThreeVector pos ;
0054     EInside       expected ;
0055 };
0056
0057 static const char* EInsideName(EInside e)
0058 {
0059     switch(e)
0060     {
0061         case kInside:  return "kInside"  ;
0062         case kSurface: return "kSurface" ;
0063         case kOutside: return "kOutside" ;
0064     }
0065     return "?" ;
0066 }
0067
0068 static int check_rays(const char* name, const G4VSolid* solid, const std::vector<Ray>& rays)
0069 {
0070     int fail = 0 ;
0071     for(const auto& r : rays)
0072     {
0073         G4double t = ssolid::Distance(solid, r.pos, r.dir, false);
0074
0075         bool miss_expected = (r.expected_t >= kInfinity) ;
0076         bool miss_got      = (t >= kInfinity) ;
0077         bool ok ;
0078
0079         if(miss_expected)
0080             ok = miss_got ;
0081         else if(miss_got)
0082             ok = false ;
0083         else
0084             ok = std::fabs(t - r.expected_t) <= r.tolerance ;
0085
0086         if(!ok)
0087         {
0088             std::cerr << "FAIL " << name << " ray=" << r.label
0089                       << "  expected=" << r.expected_t
0090                       << "  got=" << t
0091                       << "\n" ;
0092             fail++ ;
0093         }
0094         else
0095         {
0096             std::cout << "  ok " << r.label
0097                       << "  t=" << std::fixed << std::setprecision(4) << t
0098                       << "\n" ;
0099         }
0100     }
0101     return fail ;
0102 }
0103
0104 static int check_inside(const char* name, const G4VSolid* solid, const std::vector<InsidePoint>& pts)
0105 {
0106     int fail = 0 ;
0107     for(const auto& ip : pts)
0108     {
0109         EInside got = solid->Inside(ip.pos) ;
0110         bool ok = (got == ip.expected) ;
0111         if(!ok)
0112         {
0113             std::cerr << "FAIL Inside " << name << " pt=" << ip.label
0114                       << "  expected=" << EInsideName(ip.expected)
0115                       << "  got=" << EInsideName(got)
0116                       << "\n" ;
0117             fail++ ;
0118         }
0119         else
0120         {
0121             std::cout << "  ok Inside " << ip.label
0122                       << "  " << EInsideName(got)
0123                       << "\n" ;
0124         }
0125     }
0126     return fail ;
0127 }
0128
0129 // ---------------------------------------------------------------
0130 // Test 1 – OrbOrbMultiUnion1
0131 //
0132 //   3 orbs of radius 50 mm centred at x = -100, 0, +100
0133 //   (OrbOrbMultiUnion1 means num=1 → range -1..+1)
0134 //
0135 //   Geometry (Z cross-section along X axis):
0136 //
0137 //      left orb          centre orb         right orb
0138 //    x:-150..-50        x:-50..+50        x:+50..+150
0139 //
0140 //   The orbs touch at x = ±50.
0141 // ---------------------------------------------------------------
0142 static int test_OrbOrbMultiUnion()
0143 {
0144     const char* name = "OrbOrbMultiUnion1" ;
0145     const G4VSolid* solid = U4SolidMaker::Make(name);
0146     if(!solid){ std::cerr << "FAIL: could not create " << name << "\n"; return 1; }
0147
0148     const double R = 50.0 ;
0149
0150     std::vector<Ray> rays = {
0151         // Outside, shooting +Z into centre orb – should hit at z = -R → t = 200-50 = 150
0152         { "hit_centre_orb_+Z",  {0, 0, -200}, {0, 0, 1}, 150.0, kTol },
0153
0154         // Outside, shooting -Z into centre orb
0155         { "hit_centre_orb_-Z",  {0, 0,  200}, {0, 0,-1}, 150.0, kTol },
0156
0157         // Outside, shooting -X into right orb (centred at x=+100) – surface at x=+150 → t=150
0158         { "hit_right_orb_-X",   {300, 0, 0}, {-1, 0, 0}, 150.0, kTol },
0159
0160         // Outside, shooting +X into left orb (centred at x=-100) – surface at x=-150 → t=150
0161         { "hit_left_orb_+X",    {-300, 0, 0}, {1, 0, 0}, 150.0, kTol },
0162
0163         // Complete miss – ray parallel above all orbs (y > R)
0164         { "miss_above",         {0, 0, 200}, {1, 0, 0}, kInfinity, 0 },
0165
0166         // Ray from inside centre orb toward +Z surface – distance = R = 50
0167         { "inside_centre_+Z",   {0, 0, 0}, {0, 0, 1}, R, kTol },
0168
0169         // Ray from inside right orb toward +X surface – distance = R = 50
0170         { "inside_right_+X",    {100, 0, 0}, {1, 0, 0}, R, kTol },
0171
0172         // Ray along +X from far left – enters left orb at x=-150 → t = 200-150 = 50
0173         { "enter_left_along_+X", {-200, 0, 0}, {1, 0, 0}, 50.0, kTol },
0174     };
0175
0176     int fail = check_rays(name, solid, rays) ;
0177
0178     // Inside() classification checks
0179     std::vector<InsidePoint> pts = {
0180         { "origin_inside",      {0, 0, 0},       kInside  },
0181         { "right_centre",       {100, 0, 0},     kInside  },
0182         { "left_centre",        {-100, 0, 0},    kInside  },
0183         { "far_outside",        {0, 0, 500},     kOutside },
0184         { "surface_centre_+Z",  {0, 0, R},       kSurface },
0185         // x=60 is inside the right orb (centred at 100, radius 50, extends 50..150)
0186         { "inside_right_orb",   {60, 0, 0},      kInside  },
0187     };
0188
0189     fail += check_inside(name, solid, pts) ;
0190     return fail ;
0191 }
0192
0193
0194 // ---------------------------------------------------------------
0195 // Test 2 – BoxFourBoxContiguous
0196 //
0197 //   Centre box is G4Box(45, 45, 45) – a 45mm half-side cube.
0198 //   Flanking +X/-X boxes: G4Box(10, 11.5, 6.5) at x=±52
0199 //   Flanking +Y/-Y boxes: G4Box(15, 15, 6.5) at y=±50
0200 //
0201 //   Full extent: x = -62..62, y = -65..65, z = -45..45
0202 // ---------------------------------------------------------------
0203 static int test_BoxFourBoxContiguous()
0204 {
0205     const char* name = "BoxFourBoxContiguous" ;
0206     const G4VSolid* solid = U4SolidMaker::Make(name);
0207     if(!solid){ std::cerr << "FAIL: could not create " << name << "\n"; return 1; }
0208
0209     const double hz = 45.0 ; // centre box z half-extent
0210
0211     std::vector<Ray> rays = {
0212         // Shoot +Z into centre box from above → hit at z = -45 → t = 200-45 = 155
0213         { "centre_box_+Z",   {0, 0, -200}, {0, 0, 1}, 200.0 - hz, kTol },
0214
0215         // Shoot -Z into centre box from below → t = 155
0216         { "centre_box_-Z",   {0, 0,  200}, {0, 0,-1}, 200.0 - hz, kTol },
0217
0218         // Miss – ray parallel far away in Y
0219         { "miss_far_Y",      {0, 500, 0}, {0, 0, 1}, kInfinity, 0 },
0220
0221         // Ray from inside centre box → DistanceToOut(origin, +Z) = 45
0222         { "inside_centre_+Z", {0, 0, 0}, {0, 0, 1}, hz, kTol },
0223
0224         // Ray hitting +X flanking box: box extends x=42..62, z=-6.5..6.5
0225         // From (100, 0, 0) going -X → hit at x=62 → t = 38
0226         { "hit_px_flank_-X",  {100, 0, 0}, {-1, 0, 0}, 100.0 - 62.0, kTol },
0227     };
0228
0229     int fail = check_rays(name, solid, rays) ;
0230
0231     std::vector<InsidePoint> pts = {
0232         { "origin_inside",   {0, 0, 0},       kInside  },
0233         { "far_outside",     {0, 0, 500},     kOutside },
0234         { "in_px_flank",     {55, 0, 0},      kInside  },
0235     };
0236
0237     fail += check_inside(name, solid, pts) ;
0238     return fail ;
0239 }
0240
0241
0242 // ---------------------------------------------------------------
0243 // Test 3 – Verify EnsureVoxelizedMultiUnion was invoked
0244 //
0245 //   OrbOrbMultiUnion does NOT call Voxelize() itself; it relies
0246 //   on the safety net in Make().
0247 //
0248 //   Note: Geant4's voxelizer may decide to produce 0 voxels for
0249 //   a small number of solids (its internal optimisation threshold).
0250 //   So we cannot assert voxel count > 0.  Instead, we confirm the
0251 //   solid is a G4MultiUnion AND navigation returns correct results
0252 //   (which means the standard navigation path is functioning).
0253 // ---------------------------------------------------------------
0254 static int test_navigation_without_explicit_voxelize()
0255 {
0256     const char* name = "OrbOrbMultiUnion1" ;
0257     const G4VSolid* solid = U4SolidMaker::Make(name);
0258     if(!solid){ std::cerr << "FAIL: could not create " << name << "\n"; return 1; }
0259
0260     const G4MultiUnion* mu = dynamic_cast<const G4MultiUnion*>(solid) ;
0261     if(!mu){ std::cerr << "FAIL: solid is not G4MultiUnion\n"; return 1; }
0262
0263     std::cout << "  nSolids=" << mu->GetNumberOfSolids()
0264               << " voxels=" << mu->GetVoxels().GetCountOfVoxels()
0265               << "\n" ;
0266
0267     // This solid has no explicit Voxelize() call in OrbOrbMultiUnion —
0268     // EnsureVoxelizedMultiUnion in Make() should have called it.
0269     // Verify navigation correctness: from origin (inside centre orb),
0270     // DistanceToOut in +Z must equal the orb radius (50mm).
0271     G4ThreeVector pos(0, 0, 0), dir(0, 0, 1) ;
0272     EInside in ;
0273     G4double t = ssolid::Distance_(solid, pos, dir, in) ;
0274
0275     if(in != kInside)
0276     {
0277         std::cerr << "FAIL: origin should be kInside, got " << EInsideName(in) << "\n" ;
0278         return 1 ;
0279     }
0280     if(std::fabs(t - 50.0) > kTol)
0281     {
0282         std::cerr << "FAIL: expected t=50, got " << t << "\n" ;
0283         return 1 ;
0284     }
0285     std::cout << "  ok navigation correct: in=" << EInsideName(in) << " t=" << t << "\n" ;
0286     return 0 ;
0287 }
0288
0289
0290 int main(int argc, char** argv)
0291 {
0292     OPTICKS_LOG(argc, argv);
0293
0294     int fail = 0 ;
0295
0296     std::cout << "--- test_OrbOrbMultiUnion ---\n" ;
0297     fail += test_OrbOrbMultiUnion() ;
0298
0299     std::cout << "--- test_BoxFourBoxContiguous ---\n" ;
0300     fail += test_BoxFourBoxContiguous() ;
0301
0302     std::cout << "--- test_navigation_without_explicit_voxelize ---\n" ;
0303     fail += test_navigation_without_explicit_voxelize() ;
0304
0305     std::cout << "--- " << (fail == 0 ? "ALL PASS" : "FAILURES") << " ---\n" ;
0306     return fail ;
0307 }