diff --git a/CSG/CSGNode.cc b/CSG/CSGNode.cc index d3f4330a30..1a08fd9a5f 100644 --- a/CSG/CSGNode.cc +++ b/CSG/CSGNode.cc @@ -422,7 +422,10 @@ void CSGNode::setAABBLocal() { float px, py, a, radius, z1, z2 ; getParam( px, py, a, radius, z1, z2) ; - setAABB( px-radius, py-radius, z1, px+radius, py+radius, z2 ); + // F6: px,py (cx,cy) are repurposed as startPhi,deltaPhi for phi-wedge + // cylinders (sn::Cylinder 5-arg), so they are NOT the centre. Cylinder is + // origin-centred -> AABB is simply +/-radius in x/y over [z1,z2]. + setAABB(-radius, -radius, z1, radius, radius, z2); } else if( tc == CSG_DISC ) { diff --git a/CSG/csg_intersect_leaf_cylinder.h b/CSG/csg_intersect_leaf_cylinder.h index 4b9e5a3a24..f01eb7a0ec 100644 --- a/CSG/csg_intersect_leaf_cylinder.h +++ b/CSG/csg_intersect_leaf_cylinder.h @@ -33,51 +33,143 @@ into a float4 and then pick from that ? LEAF_FUNC void intersect_leaf_cylinder( bool& valid_isect, float4& isect, const quad& q0, const quad& q1, const float t_min, const float3& ray_origin, const float3& ray_direction ) { - const float& r = q0.f.w ; - const float& z1 = q1.f.x ; - const float& z2 = q1.f.y ; - const float& ox = ray_origin.x ; - const float& oy = ray_origin.y ; - const float& oz = ray_origin.z ; - const float& vx = ray_direction.x ; - const float& vy = ray_direction.y ; - const float& vz = ray_direction.z ; - - const float r2 = r*r ; - const float a = vx*vx + vy*vy ; // see CSG/sympy_cylinder.py - const float b = ox*vx + oy*vy ; - const float c = ox*ox + oy*oy - r2 ; - - float t_near, t_far, disc, sdisc ; - robust_quadratic_roots_disqualifying(t_min, t_near, t_far, disc, sdisc, a, b, c); // Solving: a t^2 + 2 b t + c = 0 - float z_near = oz+t_near*vz ; - float z_far = oz+t_far*vz ; - - const float t_z1cap = (z1 - oz)/vz ; - const float r2_z1cap = (ox+t_z1cap*vx)*(ox+t_z1cap*vx) + (oy+t_z1cap*vy)*(oy+t_z1cap*vy) ; - - const float t_z2cap = (z2 - oz)/vz ; - const float r2_z2cap = (ox+t_z2cap*vx)*(ox+t_z2cap*vx) + (oy+t_z2cap*vy)*(oy+t_z2cap*vy) ; + const float& r = q0.f.w; + const float& z1 = q1.f.x; + const float& z2 = q1.f.y; + // Optional phi-wedge clip: q0.f.x=startPhi, q0.f.y=deltaPhi (rad). + // deltaPhi==0 (or >=2π) means no clip — original full-cylinder behaviour. + const float startPhi = q0.f.x; + const float deltaPhi = q0.f.y; + const bool has_phi_clip = (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf); + const float& ox = ray_origin.x; + const float& oy = ray_origin.y; + const float& oz = ray_origin.z; + const float& vx = ray_direction.x; + const float& vy = ray_direction.y; + const float& vz = ray_direction.z; + + const float r2 = r * r; + const float a = vx * vx + vy * vy; // see CSG/sympy_cylinder.py + const float b = ox * vx + oy * vy; + const float c = ox * ox + oy * oy - r2; + + float t_near, t_far, disc, sdisc; + robust_quadratic_roots_disqualifying(t_min, t_near, t_far, disc, sdisc, a, b, c); // Solving: a t^2 + 2 b t + c = 0 + float z_near = oz + t_near * vz; + float z_far = oz + t_far * vz; + + const float t_z1cap = (z1 - oz) / vz; + const float r2_z1cap = (ox + t_z1cap * vx) * (ox + t_z1cap * vx) + (oy + t_z1cap * vy) * (oy + t_z1cap * vy); + + const float t_z2cap = (z2 - oz) / vz; + const float r2_z2cap = (ox + t_z2cap * vx) * (ox + t_z2cap * vx) + (oy + t_z2cap * vy) * (oy + t_z2cap * vy); #ifdef DEBUG - //printf("// t_z1cap %10.4f r2_z1cap %10.4f \n", t_z1cap, r2_z1cap ); - //printf("// t_z2cap %10.4f r2_z2cap %10.4f \n", t_z2cap, r2_z2cap ); + // printf("// t_z1cap %10.4f r2_z1cap %10.4f \n", t_z1cap, r2_z1cap ); + // printf("// t_z2cap %10.4f r2_z2cap %10.4f \n", t_z2cap, r2_z2cap ); #endif - float t_cand = CUDART_INF_F ; - if( t_near > t_min && z_near > z1 && z_near < z2 && t_near < t_cand ) t_cand = t_near ; - if( t_far > t_min && z_far > z1 && z_far < z2 && t_far < t_cand ) t_cand = t_far ; - if( t_z1cap > t_min && r2_z1cap <= r2 && t_z1cap < t_cand ) t_cand = t_z1cap ; - if( t_z2cap > t_min && r2_z2cap <= r2 && t_z2cap < t_cand ) t_cand = t_z2cap ; - - valid_isect = t_cand > t_min && t_cand < CUDART_INF_F ; + auto phi_in_wedge = [&](float xh, float yh) -> bool { + if (!has_phi_clip) + return true; + float phi = atan2f(yh, xh); + const float twoPi = 2.f * M_PIf; + while (phi < startPhi) phi += twoPi; + return (phi - startPhi) <= deltaPhi; + }; + // Side-surface candidates also need to be inside the wedge. + bool t_near_ok = (z_near > z1 && z_near < z2) && phi_in_wedge(ox + t_near * vx, oy + t_near * vy); + bool t_far_ok = (z_far > z1 && z_far < z2) && phi_in_wedge(ox + t_far * vx, oy + t_far * vy); + // Cap candidates need their hit point's r AND phi inside. + bool t_z1cap_ok = (r2_z1cap <= r2) && phi_in_wedge(ox + t_z1cap * vx, oy + t_z1cap * vy); + bool t_z2cap_ok = (r2_z2cap <= r2) && phi_in_wedge(ox + t_z2cap * vx, oy + t_z2cap * vy); + + // Flat radial wedge-wall candidates at phi=startPhi and phi=startPhi+deltaPhi. + // The two half-planes contain the z-axis. A hit must land within + // r_h^2 <= r^2, z1 < z_h < z2, and on the +r side of the wall + // (i.e. not the antipodal half-plane through the same line). + // Outward normal at startPhi: n1 = ( sin(startPhi), -cos(startPhi), 0) + // Outward normal at startPhi+deltaPhi: n2 = (-sin(endPhi), cos(endPhi), 0) + float t_phi1 = CUDART_INF_F, t_phi2 = CUDART_INF_F; + bool t_phi1_ok = false, t_phi2_ok = false; + float n1x = 0.f, n1y = 0.f, n2x = 0.f, n2y = 0.f; +#ifndef CSG_CYL_NO_PHI_WALLS + if (has_phi_clip) + { + const float endPhi = startPhi + deltaPhi; + const float s1 = sinf(startPhi), c1 = cosf(startPhi); + const float s2 = sinf(endPhi), c2 = cosf(endPhi); + n1x = s1; + n1y = -c1; + n2x = -s2; + n2y = c2; + + // phi=startPhi plane: { x*s1 - y*c1 = 0 } solved for t + const float denom1 = vx * s1 - vy * c1; + if (fabsf(denom1) > 1.e-12f) + { + t_phi1 = -(ox * s1 - oy * c1) / denom1; + const float xh = ox + t_phi1 * vx; + const float yh = oy + t_phi1 * vy; + const float zh = oz + t_phi1 * vz; + const float r2_h = xh * xh + yh * yh; + const bool on_pos_r = (xh * c1 + yh * s1) >= 0.f; // not antipodal half + t_phi1_ok = on_pos_r && (r2_h <= r2) && (zh > z1) && (zh < z2); + } + // phi=endPhi plane: { x*s2 - y*c2 = 0 } solved for t + const float denom2 = vx * s2 - vy * c2; + if (fabsf(denom2) > 1.e-12f) + { + t_phi2 = -(ox * s2 - oy * c2) / denom2; + const float xh = ox + t_phi2 * vx; + const float yh = oy + t_phi2 * vy; + const float zh = oz + t_phi2 * vz; + const float r2_h = xh * xh + yh * yh; + const bool on_pos_r = (xh * c2 + yh * s2) >= 0.f; + t_phi2_ok = on_pos_r && (r2_h <= r2) && (zh > z1) && (zh < z2); + } + } +#endif // CSG_CYL_NO_PHI_WALLS + + float t_cand = CUDART_INF_F; + if (t_near > t_min && t_near_ok && t_near < t_cand) + t_cand = t_near; + if (t_far > t_min && t_far_ok && t_far < t_cand) + t_cand = t_far; + if (t_z1cap > t_min && t_z1cap_ok && t_z1cap < t_cand) + t_cand = t_z1cap; + if (t_z2cap > t_min && t_z2cap_ok && t_z2cap < t_cand) + t_cand = t_z2cap; + if (t_phi1 > t_min && t_phi1_ok && t_phi1 < t_cand) + t_cand = t_phi1; + if (t_phi2 > t_min && t_phi2_ok && t_phi2 < t_cand) + t_cand = t_phi2; + + valid_isect = t_cand > t_min && t_cand < CUDART_INF_F; if(valid_isect) { - bool sheet = ( t_cand == t_near || t_cand == t_far ) ; - isect.x = sheet ? (ox + t_cand*vx)/r : 0.f ; - isect.y = sheet ? (oy + t_cand*vy)/r : 0.f ; - isect.z = sheet ? 0.f : ( t_cand == t_z1cap ? -1.f : 1.f) ; - isect.w = t_cand ; + bool sheet = (t_cand == t_near || t_cand == t_far); + bool cap = (t_cand == t_z1cap || t_cand == t_z2cap); + bool phi1 = (t_cand == t_phi1); + if (sheet) + { + isect.x = (ox + t_cand * vx) / r; + isect.y = (oy + t_cand * vy) / r; + isect.z = 0.f; + } + else if (cap) + { + isect.x = 0.f; + isect.y = 0.f; + isect.z = (t_cand == t_z1cap ? -1.f : 1.f); + } + else // phi-wall + { + isect.x = phi1 ? n1x : n2x; + isect.y = phi1 ? n1y : n2y; + isect.z = 0.f; + } + isect.w = t_cand; } } @@ -136,18 +228,34 @@ The SDF rules for CSG combinations:: LEAF_FUNC float distance_leaf_cylinder( const float3& pos, const quad& q0, const quad& q1 ) { - const float radius = q0.f.w ; - const float z1 = q1.f.x ; - const float z2 = q1.f.y ; // z2 > z1 + const float radius = q0.f.w; + const float z1 = q1.f.x; + const float z2 = q1.f.y; // z2 > z1 + const float startPhi = q0.f.x; + const float deltaPhi = q0.f.y; - float sd_capslab = fmaxf( pos.z - z2 , z1 - pos.z ); - float sd_infcyl = sqrtf( pos.x*pos.x + pos.y*pos.y ) - radius ; - float sd = fmaxf( sd_capslab, sd_infcyl ); + float sd_capslab = fmaxf(pos.z - z2, z1 - pos.z); + float sd_infcyl = sqrtf(pos.x * pos.x + pos.y * pos.y) - radius; + float sd = fmaxf(sd_capslab, sd_infcyl); + + if (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf) + { + float phi = atan2f(pos.y, pos.x); + const float twoPi = 2.f * M_PIf; + while (phi < startPhi) phi += twoPi; + float dphi = phi - startPhi; + if (dphi > deltaPhi) + { + float r_xy = sqrtf(pos.x * pos.x + pos.y * pos.y); + float d_to_edge = r_xy * fminf(dphi - deltaPhi, twoPi - dphi); + sd = fmaxf(sd, d_to_edge); + } + } #ifdef DEBUG - printf("//distance_leaf_cylinder sd %10.4f \n", sd ); + printf("//distance_leaf_cylinder sd %10.4f \n", sd); #endif - return sd ; + return sd; } diff --git a/CSG/csg_intersect_leaf_zsphere.h b/CSG/csg_intersect_leaf_zsphere.h index ef7d6b0a62..eeb2e200d3 100644 --- a/CSG/csg_intersect_leaf_zsphere.h +++ b/CSG/csg_intersect_leaf_zsphere.h @@ -3,18 +3,38 @@ LEAF_FUNC float distance_leaf_zsphere(const float3& pos, const quad& q0, const quad& q1 ) { - float3 center = make_float3(q0.f); - float radius = q0.f.w; + // ZSphere is always centred at the local origin (G4Sphere doesn't carry + // a translation in its primitive — that's done by the parent transform). + // q0.f.x and q0.f.y therefore carry an OPTIONAL phi-wedge clip: + // q0.f.x = startPhi (rad), q0.f.y = deltaPhi (rad). + // deltaPhi == 0 (or ≥ 2π) means no clip — original ZSphere behaviour. + const float startPhi = q0.f.x; + const float deltaPhi = q0.f.y; + const float radius = q0.f.w; const float2 zdelta = make_float2(q1.f); - const float z2 = center.z + zdelta.y ; - const float z1 = center.z + zdelta.x ; + const float z2 = zdelta.y; + const float z1 = zdelta.x; - float3 p = pos - center; - float sd_sphere = length(p) - radius ; - float sd_capslab = fmaxf( pos.z - z2 , z1 - pos.z ); + const float3 p = pos; + float sd_sphere = length(p) - radius; + float sd_capslab = fmaxf(pos.z - z2, z1 - pos.z); float sd = fmaxf( sd_capslab, sd_sphere ); // CSG intersect - return sd ; + + if (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf) + { + float phi = atan2f(p.y, p.x); + const float twoPi = 2.f * M_PIf; + while (phi < startPhi) phi += twoPi; + float dphi = phi - startPhi; + if (dphi > deltaPhi) + { + float r_xy = sqrtf(p.x * p.x + p.y * p.y); + float d_to_edge = r_xy * fminf(dphi - deltaPhi, twoPi - dphi); + sd = fmaxf(sd, d_to_edge); + } + } + return sd; } /** @@ -44,10 +64,24 @@ See : notes/issues/unexpected_zsphere_miss_from_inside_for_rays_that_would_be_ex LEAF_FUNC void intersect_leaf_zsphere(bool& valid_isect, float4& isect, const quad& q0, const quad& q1, const float& t_min, const float3& ray_origin, const float3& ray_direction ) { - const float3 center = make_float3(q0.f); - float3 O = ray_origin - center; + // ZSphere is always centred at the local origin; q0.f.x/q0.f.y carry an + // OPTIONAL phi-wedge clip (see distance_leaf_zsphere comment). + const float3 center = make_float3(0.f, 0.f, 0.f); + float3 O = ray_origin - center; float3 D = ray_direction; const float radius = q0.f.w; +#ifdef ZSPHERE_PHI_DEBUG + if (q0.f.y > 0.f && q0.f.y < 6.28f) + { + static int cnt = 0; + if (cnt < 5) + { + printf("// ZSphere phi: startPhi=%.4f deltaPhi=%.4f radius=%.2f z1=%.2f z2=%.2f\n", + q0.f.x, q0.f.y, q0.f.w, q1.f.x, q1.f.y); + cnt++; + } + } +#endif float b = dot(O, D); // t of closest approach to sphere center float c = dot(O, O)-radius*radius; // < 0. indicates ray_origin inside sphere @@ -109,18 +143,41 @@ void intersect_leaf_zsphere(bool& valid_isect, float4& isect, const quad& q0, co // hmm somehow is seems unclean to have to use both z and t language - float t_cand = t_min ; + // Optional phi-wedge clip on the sphere-surface candidates. + // q0.f.x = startPhi (rad), q0.f.y = deltaPhi (rad). deltaPhi == 0 + // (or ≥ 2π) means no clip — the original two-cap ZSphere behaviour. + const float startPhi = q0.f.x; + const float deltaPhi = q0.f.y; + const bool has_phi_clip = (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf); + auto phi_in_wedge = [&](float t) -> bool { + if (!has_phi_clip) + return true; + float x = O.x + t * D.x; // local-frame hit coords (O = origin − center) + float y = O.y + t * D.y; + float phi = atan2f(y, x); + const float twoPi = 2.f * M_PIf; + while (phi < startPhi) phi += twoPi; + float dphi = phi - startPhi; + return dphi <= deltaPhi; + }; + bool t1sph_ok = (z1sph > zmin && z1sph <= zmax) && phi_in_wedge(t1sph); + bool t2sph_ok = (z2sph > zmin && z2sph <= zmax) && phi_in_wedge(t2sph); + + float t_cand = t_min; if(sdisc > 0.f) { #ifdef DEBUG_RECORD - //std::raise(SIGINT); + // std::raise(SIGINT); #endif - if( t1sph > t_min && z1sph > zmin && z1sph <= zmax ) t_cand = t1sph ; // t1sph qualified and t1cap disabled or disqualified -> t1sph - else if( t1cap > t_min ) t_cand = t1cap ; // t1cap qualifies -> t1cap + if (t1sph > t_min && t1sph_ok) + t_cand = t1sph; // t1sph qualified and t1cap disabled or disqualified -> t1sph + else if (t1cap > t_min) + t_cand = t1cap; // t1cap qualifies -> t1cap else if( t2cap > t_min ) t_cand = t2cap ; // t2cap qualifies -> t2cap - else if( t2sph > t_min && z2sph > zmin && z2sph <= zmax) t_cand = t2sph ; // t2sph qualifies and t2cap disabled or disqialified -> t2sph + else if (t2sph > t_min && t2sph_ok) + t_cand = t2sph; // t2sph qualifies and t2cap disabled or disqialified -> t2sph } valid_isect = t_cand > t_min ; diff --git a/CSG/tests/CMakeLists.txt b/CSG/tests/CMakeLists.txt index fa0692d703..ea5cb79f8d 100644 --- a/CSG/tests/CMakeLists.txt +++ b/CSG/tests/CMakeLists.txt @@ -43,7 +43,8 @@ set(TEST_SOURCES CSGCopyTest.cc - intersect_leaf_phicut_test.cc + intersect_leaf_phicut_test.cc + intersect_leaf_phi_wedge_test.cc intersect_leaf_thetacut_test.cc intersect_leaf_box3_test.cc intersect_leaf_cylinder_test.cc diff --git a/CSG/tests/intersect_leaf_phi_wedge_test.cc b/CSG/tests/intersect_leaf_phi_wedge_test.cc new file mode 100644 index 0000000000..6c412c1044 --- /dev/null +++ b/CSG/tests/intersect_leaf_phi_wedge_test.cc @@ -0,0 +1,164 @@ +/** +intersect_leaf_phi_wedge_test.cc +================================== + +Self-checking regression test for the phi-wedge clip baked into the ZSphere +(csg_intersect_leaf_zsphere.h) and Cylinder (csg_intersect_leaf_cylinder.h) +leaf primitives. + +Unlike intersect_leaf_phicut_test.cc (which exercises the standalone phicut +half-space node), this targets the per-primitive phi-wedge encoding: + + q0.f.x = startPhi (rad) q0.f.y = deltaPhi (rad) q0.f.w = radius + q1.f.x = z1 q1.f.y = z2 + +It pins the behaviour the original CSG_INTERSECTION(prim, PhiCut) leak got +wrong: a ray whose NEAR surface candidate falls outside the wedge must still +find the in-wedge surface/wall further along the ray (rather than be lost), +while a ray that crosses the primitive entirely outside the wedge must MISS. + +CPU-only: the leaf intersect headers are shared host/device code. +Exit code is non-zero if any assertion fails, so it is meaningful under ctest. +**/ + +#include +#include + +#include "scuda.h" +#include "sqat4.h" +#include "squad.h" + +// required include order: csg_intersect_leaf_head.h defines LEAF_FUNC etc. +// used by csg_robust_quadratic_roots.h and the leaf intersect headers below +// clang-format off +#include "csg_intersect_leaf_head.h" +#include "csg_robust_quadratic_roots.h" +#include "csg_intersect_leaf_zsphere.h" +#include "csg_intersect_leaf_cylinder.h" +// clang-format on + +static int g_fail = 0; + +static void check(const char* tag, bool cond) +{ + printf("// %-58s %s\n", tag, cond ? "PASS" : "*** FAIL ***"); + if (!cond) + g_fail += 1; +} + +// relative-ish float compare, generous enough for fp32 geometry +static bool feq(float a, float b, float eps = 2e-2f) +{ + return fabsf(a - b) <= eps * (1.f + fabsf(b)); +} + +static quad q_phi(float startPhi, float deltaPhi, float radius) +{ + quad q; + q.f = make_float4(startPhi, deltaPhi, 0.f, radius); + return q; +} +static quad q_zrange(float z1, float z2) +{ + quad q; + q.f = make_float4(z1, z2, 0.f, 0.f); + return q; +} + +static bool run_zsphere(const quad& q0, const quad& q1, float3 o, float3 d, float3& hit, float& t) +{ + float4 isect = make_float4(0.f, 0.f, 0.f, 0.f); + bool valid = false; + intersect_leaf_zsphere(valid, isect, q0, q1, 0.f, o, d); + t = isect.w; + hit = make_float3(o.x + t * d.x, o.y + t * d.y, o.z + t * d.z); + return valid; +} +static bool run_cylinder(const quad& q0, const quad& q1, float3 o, float3 d, float3& hit, float& t) +{ + float4 isect = make_float4(0.f, 0.f, 0.f, 0.f); + bool valid = false; + intersect_leaf_cylinder(valid, isect, q0, q1, 0.f, o, d); + t = isect.w; + hit = make_float3(o.x + t * d.x, o.y + t * d.y, o.z + t * d.z); + return valid; +} + +// azimuth of a hit point, normalized into [startPhi, startPhi+2pi) +static bool phi_in_wedge(float3 h, float startPhi, float deltaPhi) +{ + float phi = atan2f(h.y, h.x); + while (phi < startPhi) + phi += 2.f * M_PIf; + return (phi - startPhi) <= deltaPhi + 1e-4f; +} + +void test_zsphere_phi() +{ + const float R = 100.f; + const float HALF = 0.5f * M_PIf; // deltaPhi = pi/2 wedge [0, pi/2] + const quad q1 = q_zrange(-R, R); // full sphere in z + const quad qw = q_phi(0.f, HALF, R); + float3 h; + float t; + const float c = cosf(0.25f * M_PIf), s = sinf(0.25f * M_PIf); + + // (1) ray aimed at a surface point squarely inside the wedge (phi=pi/4) -> hit + bool v1 = run_zsphere(qw, q1, make_float3(2 * R * c, 2 * R * s, 0.f), make_float3(-c, -s, 0.f), h, t); + check("zsphere: in-wedge ray hits", v1); + check("zsphere: in-wedge hit on surface inside wedge", v1 && feq(h.x, R * c) && feq(h.y, R * s) && phi_in_wedge(h, 0.f, HALF)); + + // (2) LEAK regression: near surface candidate is OUTSIDE the wedge (phi~150deg), + // the far candidate (phi~30deg) is inside -> must be recovered, not lost. + bool v2 = run_zsphere(qw, q1, make_float3(-2 * R, 0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("zsphere: leak-case recovers far in-wedge hit (was lost)", v2); + check("zsphere: leak-case hit is the +x (far) surface, in wedge", v2 && h.x > 0.f && feq(h.y, 0.5f * R) && phi_in_wedge(h, 0.f, HALF)); + + // (3) ray crossing the sphere entirely outside the wedge (y<0 side) -> MISS + bool v3 = run_zsphere(qw, q1, make_float3(-2 * R, -0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("zsphere: ray fully outside wedge misses", !v3); + + // (4) control: same geometry as (2) but NO phi clip (deltaPhi=0) -> takes the NEAR (-x) surface + bool v4 = run_zsphere(q_phi(0.f, 0.f, R), q1, make_float3(-2 * R, 0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("zsphere: no-clip control hits", v4); + check("zsphere: no-clip control takes near (-x) surface", v4 && h.x < 0.f && feq(h.y, 0.5f * R)); +} + +void test_cylinder_phi() +{ + const float R = 100.f; + const float HALF = 0.5f * M_PIf; // [0, pi/2] + const quad q1 = q_zrange(-R, R); + const quad qw = q_phi(0.f, HALF, R); + float3 h; + float t; + const float c = cosf(0.25f * M_PIf), s = sinf(0.25f * M_PIf); + + // (1) ray aimed at the curved surface inside the wedge (phi=pi/4) -> hit + bool v1 = run_cylinder(qw, q1, make_float3(2 * R * c, 2 * R * s, 0.f), make_float3(-c, -s, 0.f), h, t); + check("cylinder: in-wedge ray hits curved surface", v1); + check("cylinder: in-wedge hit inside wedge & on radius", v1 && feq(h.x, R * c) && feq(h.y, R * s) && phi_in_wedge(h, 0.f, HALF)); + + // (2) wedge-boundary regression: ray enters from -x at y=+50; the near curved + // candidate (phi~150deg) is clipped, so the first valid hit must be in-wedge + // (the phi=pi/2 wall at x=0, or the far curved surface) -- NOT a leak-through. + bool v2 = run_cylinder(qw, q1, make_float3(-2 * R, 0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("cylinder: wedge-boundary ray yields an in-wedge hit (no leak)", v2); + check("cylinder: wedge-boundary hit is inside the wedge", v2 && phi_in_wedge(h, 0.f, HALF) && h.y > 0.f); + + // (3) ray crossing entirely outside the wedge (y<0) -> MISS + bool v3 = run_cylinder(qw, q1, make_float3(-2 * R, -0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("cylinder: ray fully outside wedge misses", !v3); + + // (4) control: NO phi clip -> near (-x) curved surface + bool v4 = run_cylinder(q_phi(0.f, 0.f, R), q1, make_float3(-2 * R, 0.5f * R, 0.f), make_float3(1.f, 0.f, 0.f), h, t); + check("cylinder: no-clip control hits near (-x) surface", v4 && h.x < 0.f && feq(h.y, 0.5f * R)); +} + +int main(int argc, char** argv) +{ + test_zsphere_phi(); + test_cylinder_phi(); + printf("// intersect_leaf_phi_wedge_test : %s (%d failure%s)\n", g_fail == 0 ? "PASS" : "FAIL", g_fail, g_fail == 1 ? "" : "s"); + return g_fail == 0 ? 0 : 1; +} diff --git a/sysrap/sn.h b/sysrap/sn.h index 327400d59b..c0cdccbb3e 100644 --- a/sysrap/sn.h +++ b/sysrap/sn.h @@ -441,6 +441,7 @@ struct SYSRAP_API sn std::string descXF() const ; static sn* Cylinder(double radius, double z1, double z2) ; + static sn* Cylinder(double radius, double z1, double z2, double startPhi, double deltaPhi); // F5 phi-baked static sn* CutCylinder( double R, double dz, @@ -468,6 +469,7 @@ struct SYSRAP_API sn static sn* Cone(double r1, double z1, double r2, double z2); static sn* Sphere(double radius); static sn* ZSphere(double radius, double z1, double z2); + static sn* ZSphere(double radius, double z1, double z2, double startPhi, double deltaPhi); // F4 phi-baked static sn* Box3(double fullside); static sn* Box3(double fx, double fy, double fz ); static sn* Torus(double rmin, double rmax, double rtor, double startPhi_deg, double deltaPhi_deg ); @@ -3062,6 +3064,17 @@ inline sn* sn::Cylinder(double radius, double z1, double z2) // static nd->setBB( -radius, -radius, z1, +radius, +radius, z2 ); return nd ; } +inline sn* sn::Cylinder(double radius, double z1, double z2, double startPhi, double deltaPhi) // static +{ + assert(z2 > z1); + sn* nd = Create(CSG_CYLINDER); + // F5 phi-bake: park startPhi/deltaPhi (radians) in cx/cy slots (0 for a + // centred cylinder). csg_intersect_leaf_cylinder.h treats deltaPhi==0 + // (or >=2pi) as "no clip". Same convention as the phi-aware ZSphere. + nd->setPA(startPhi, deltaPhi, 0.f, radius, z1, z2); + nd->setBB(-radius, -radius, z1, +radius, +radius, z2); + return nd; +} /** sn::CutCylinder @@ -3242,6 +3255,19 @@ inline sn* sn::ZSphere(double radius, double z1, double z2) // static nd->setBB( -radius, -radius, z1, radius, radius, z2 ); return nd ; } +inline sn* sn::ZSphere(double radius, double z1, double z2, double startPhi, double deltaPhi) // static +{ + assert(radius > zero); + assert(z2 > z1); + sn* nd = Create(CSG_ZSPHERE); + // F4 phi-bake: park startPhi/deltaPhi (radians) into the cx/cy param slots + // (always zero for a centred ZSphere, so safely repurposed). q1.f.z/.w are + // reserved by CSGNode for boundary/index. csg_intersect_leaf_zsphere.h treats + // deltaPhi==0 (or >=2pi) as "no clip", so the 3-arg ZSphere is unchanged. + nd->setPA(startPhi, deltaPhi, zero, radius, z1, z2); + nd->setBB(-radius, -radius, z1, radius, radius, z2); + return nd; +} inline sn* sn::Box3(double fullside) // static { return Box3(fullside, fullside, fullside); diff --git a/u4/U4Solid.h b/u4/U4Solid.h index 9a47619d58..0edccabedd 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -556,11 +556,21 @@ inline sn* U4Solid::init_Sphere_(char layer) double deltaPhi = sphere->GetDeltaPhiAngle()/CLHEP::radian ; bool has_deltaPhi = startPhi != 0. || deltaPhi != 2.*CLHEP::pi ; - bool has_deltaPhi_expect = has_deltaPhi == false ; - assert( has_deltaPhi_expect ); - if(!has_deltaPhi_expect) std::raise(SIGINT); - - return z_slice ? sn::ZSphere( radius, zmin, zmax ) : sn::Sphere(radius ) ; + // F4 phi-bake: a phi-wedge G4Sphere now builds a phi-clipped ZSphere + // (startPhi/deltaPhi baked into the primitive) instead of aborting. + // See sn::ZSphere(5-arg) + csg_intersect_leaf_zsphere.h. + if (z_slice || has_deltaPhi) + { + if (!z_slice) + { + zmin = -radius; + zmax = +radius; + } + if (has_deltaPhi) + return sn::ZSphere(radius, zmin, zmax, startPhi, deltaPhi); + return sn::ZSphere(radius, zmin, zmax); + } + return sn::Sphere(radius); } @@ -800,8 +810,19 @@ inline void U4Solid::init_Tubs() double rmin = tubs->GetInnerRadius()/CLHEP::mm ; bool has_inner = rmin > 0. ; - sn* outer = sn::Cylinder(rmax, -hz, hz ); + // F5 phi-bake: build phi-clipped cylinders for a G4Tubs with a phi wedge + // (startPhi/deltaPhi baked into the primitive via sn::Cylinder 5-arg). + double startPhi = tubs->GetStartPhiAngle() / CLHEP::radian; + double deltaPhi = tubs->GetDeltaPhiAngle() / CLHEP::radian; + bool has_deltaPhi = startPhi != 0. || deltaPhi != 2. * CLHEP::pi; + + auto makeCyl = [&](double R, double zlo, double zhi) { + return has_deltaPhi + ? sn::Cylinder(R, zlo, zhi, startPhi, deltaPhi) + : sn::Cylinder(R, zlo, zhi); + }; + sn* outer = makeCyl(rmax, -hz, hz); if(has_inner == false) { @@ -813,7 +834,7 @@ inline void U4Solid::init_Tubs() double nudge_inner = 0.01 ; double dz = do_nudge_inner ? hz*nudge_inner : 0. ; - sn* inner = sn::Cylinder(rmin, -(hz+dz), hz+dz ); + sn* inner = makeCyl(rmin, -(hz + dz), hz + dz); root = sn::Boolean( CSG_DIFFERENCE, outer, inner ); }