Created attachment 492 [details] Minimal test case with plots The attached minimal example demonstrates a problem with Geant4's calculations of normal vectors to the surface of cones. As a reference, I am testing the same geometry with ROOT. Compile the test file with g++ ../geomTest.cc `geant4-config --cflags --libs` `root-config --cflags --evelibs` -o geomTest The test is the calculation of the surface normal around the following point (units in cm): const G4double x = -1.331514e+00, y=-9.139028e-01, z=-2.876654e-01; The test program yields (amended for readability): Geant4 TEST dist = 5.85527e-07 normal = (0.506501,0.493413,-0.707107) ROOT TEST dist = 5.85527e-07 normal = (0.858151,0.141849,0.493413) Clearly Geant4 and ROOT agree on the distance from the surface, but not on the direction of the normal vector. Note that the point is not close to any of the corners of the shape. I believe ROOT's answer is the correct one. I am attaching 2D sections of the geometry around the considered point, which show that all the components of the sought normal (which should point OUT from "TOP" and INTO "cone") should be positive. Tested with Geant4 v10.4.
We cannot reproduce the problem. The test below int main() { const double coneRadius = 4.; const double coneHeight = 4.; G4Cons *solid = new G4Cons("cone", 0.*cm, coneRadius*cm, // rMin1, rMax1 0.*cm, 0.*cm, // rMin2, rMax2 0.5*coneHeight*cm, // dZ 0., twopi); // phiMin, deltaPhi G4ThreeVector p(-1.331514e+00*cm, -9.139028e-01*cm, -2.876654e-01*cm); G4ThreeVector v(4.215310e-01, 8.672989e-01, 2.647721e-01); G4double distToIn = solid->DistanceToIn(p); G4double distToOut = solid->DistanceToOut(p); G4cout << "distToIn, distToOut = " << distToIn << ", " << distToOut << G4endl; G4ThreeVector normal = solid->SurfaceNormal(p); G4cout << "normal = " << normal << G4endl; G4bool trueNorm; double dist = solid->DistanceToOut(p,v, true, &trueNorm, &normal); G4cout << "dist = " << dist << ", normal = " << normal << G4endl; } produces the following output distToIn, distToOut = 0, 4.75663 normal = (-0.582995,-0.400146,0.707107) dist = 28.2228, normal = (-0.0651076,0.704103,0.707107) The lateral surface of the cone is inclined at 45 degrees to the z-plane, so the z component of the normal should be equal to 0.707107. That perfectly corresponds to the output above.
Evgueni, thank you for looking into this. I think the rotation part (which you omitted from your test) is crucial: the bug sometimes disappears for certain values of the Euler angles for the cone rotation. Some triples of (phi, theta, psi) angles that trigger the bug here: (-45, 90, 0) (0, 90, 0) // x component is correct, y and z are off by a 90° turn (-30, 90, 0) Some triples of (phi, theta, psi) angles that *do not* trigger the bug here: (0, 0, 0) (30, 90, 0) (30, 90, 30) (30, 30, 30) Also, I can trigger the bug if I replace the cone with a cylinder (G4Tubs), but (as far as I can tell) not with boxes (G4Box) or spheres (G4Sphere).
Dear Davide, It looks that the source of the problem is a misuse of navigator->GetGlobalExitNormal(point, &valid), and not a bug in G4Cons (or G4Tubs). The comment to the G4Navigator::GetGlobalExitNormal() says: // Return Exit Surface Normal and validity too. // Can only be called if the Navigator's last Step has crossed a // volume geometrical boundary. // It returns the Normal to the surface pointing out of the volume that // was left behind and/or into the volume that was entered. In you case the point did not cross any boundary. After navigator->LocateGlobalPointAndSetup(point) the points is in the TOP volume. Then dist = navigator->ComputeStep(point, omega, kInfinity, safety) returns the distance to the cone. And just after that, there is a call normal = navigator->GetGlobalExitNormal(point, &valid) I've included some prints out into the function to understand what happened inside it. The function has calculated a local normal which points INSIDE the cone. The normal looks ok, z-value is equal to -0.707107 (inclination in 45 degrees). But then, an IDENTITY transformation has been applied to convert the local normal to the global normal. It means that a conversion matrix has not been set: ==== Enter GetGlobalExitNormal ==== localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal) ==== globalNormal = localToGlobal.TransformAxis( localNormal ) ==== localNormal = (0.506501,0.493413,-0.707107) ==== globalNormal = (0.506501,0.493413,-0.707107) I added to the conversation John Apostolakis. He understand better then me how the function should be used, and possibly can advice you on how to fix the test. Best regards, Evgueni
I've forgotten to mention that if to invert the rotation matrix defined in your test and apply it to the normal: (0.506501,0.493413,-0.707107) then the result will be as expected: (0.85815,0.14185,0.493413)
Evgueni, I'm confused. I thought the call to ComputeStep() was sufficient to satisfy the preconditions to GetGlobalExitNormal(). I had noticed the same thing about the rotation matrix. Let's see if John can shed some light on this matter. Thanks again, Davide
One more detail: there may or may not be a small bug in the example code that I provided, insofar as GetGlobalExitNormal() is called on the starting point. I replaced the call with G4ThreeVector normal = navigator->GetGlobalExitNormal(point + dist * omega, &valid); However, this results in the same normal.
Hi Davide, I think the observation Evgueni has made is correct. ComputeStep() does NOT move/relocate your point to the new position, so the call you do to GetGlobalExitNormal() is not valid (see the explanation in the comments for that method in G4Navigator). The function can only be called if the Navigator's last Step has crossed the cone's boundary. But John can certainly better explain this. Gabriele
Hi Gabriele, So should I call LocateGlobalPointAndSetup() on the new point after ComputeStep()? What is the correct way to satisfy the precondition on GetGlobalExitNormal()? Thank you, Davide
Me again, I looked at G4Navigator.cc and it turns out that ComputeStep() *does* update fStepEndPoint (around line 1060): fStepEndPoint = pGlobalpoint + std::min(Step,pCurrentProposedStepLength) * pDirection; Also, I tried adding a second call to LocateGlobalPointAndSetup() or a call to LocateGlobalPointWithinVolume() before GetGlobalExitNormal(), but Geant4 did not like it: -------- WWWW ------- G4Exception-START -------- WWWW ------- *** ExceptionHandler is not defined *** *** G4Exception : GeomNav0003 issued by : G4Navigator::GetLocalExitNormal() Function called when *NOT* at a Boundary. Exit Normal not calculated. *** This is just a warning message. *** -------- WWWW -------- G4Exception-END --------- WWWW ------- Geant4 TEST dist = 5.85527e-07 normal = (0,0,0) Note that no such warning is emitted by the original test code, which seems to suggest that the preconditions are satisfied after all. So I'm stumped again.
I confirm what Davide has said - when G4Navigator::ComputeStep is called, the endpoint of that step is prepared for a call to GetGlobalExitNormal(). This functionality is already used in the propagation of field, so it should work!
Thanks John for clarifying this. So, sorry, no issues with the pre-conditions for calling GetGlobalExitNormal() from G4Navigator. Looking at the specifications for GetGlobalExitNormal() in G4Navigator.hh, one reads: // Return Exit Surface Normal and validity too. // Can only be called if the Navigator's last Step has crossed a // volume geometrical boundary. // It returns the Normal to the surface pointing out of the volume that // was left behind and/or into the volume that was entered. Therefore, the normal returned is in -local- coordinates system (i.e. no rotation is applied in this case, and, since the particle is entering the cone, the normal points inside the solid). According to what also reported by Evgueni, seems to me everything is behaving as expected (i.e. the value returned by ROOT in global coordinates can be obtained by applying the inverse rotation)! I don't know what are the conventions used in ROOT, but I believe the test is comparing different things..
Sorry, I must take this back. GetGlobalExitNormal() should return the normal vector in -global- coordinates. If this is not the case and the rotation-matrix is not
Sorry, I must take this back. GetGlobalExitNormal() should return the normal vector in -global- coordinates. If this is not the case and the rotation-matrix is not taken into consideration, then there is a real problem. Reopening the ticket and assigning it to navigation.
I have modified the test to correctly take into account units and call GetGlobalExitNormal() on the new point: void g4test(G4double const x, G4double const y, G4double const z, G4double const vx, G4double const vy, G4double const vz) { G4Navigator *navigator = new G4Navigator; G4VPhysicalVolume *world = g4geometry(); navigator->SetWorldVolume(world); G4ThreeVector point(x, y, z); G4ThreeVector omega(vx, vy, vz); G4VPhysicalVolume* vol = navigator->LocateGlobalPointAndSetup(point); G4double safety = 0.0; G4double dist = navigator->ComputeStep(point, omega, kInfinity, safety); G4bool valid = false; G4ThreeVector new_point = point + dist*omega; // G4VPhysicalVolume* new_vol = navigator->LocateGlobalPointAndSetup(new_point); G4ThreeVector normal = navigator->GetGlobalExitNormal(new_point, &valid); G4cout.precision(16); G4cout << "Geant4 TEST\n" << " first volume = " << vol->GetName() << '\n' // << " last volume = " << new_vol->GetName() << '\n' << " point = " << point << '\n' << " new point = " << new_point << '\n' << " dist = " << dist/cm << '\n' << " normal = " << normal << '\n' << " valid = " << valid << G4endl; } int main() { const G4double x = -1.331514e+00*cm, y=-9.139028e-01*cm, z=-2.876654e-01*cm; const G4double vx = 4.215310e-01, vy=8.672989e-01, vz=2.647721e-01; g4test(x, y, z, vx, vy, vz); return 0; } I still get the wrong result (and the point is still located on the Box): Checking overlaps for volume cone ... OK! Geant4 TEST first volume = TOP point = (-13.31514,-9.139028,-2.876654) new point = (-13.31513753182348,-9.139022921733442,-2.876652449688682) dist = 5.855266917542679e-07 normal = (0.5065013452049695,0.4934129987196894,-0.7071067811865475) valid = 1 If I uncomment the call to LocateGlobalPointAndSetup() on the new point, then I get the expected result (and the point is located on the Cone)! Checking overlaps for volume cone ... OK! Geant4 TEST first volume = TOP last volume = cone point = (-13.31514,-9.139028,-2.876654) new point = (-13.31513753182348,-9.139022921733442,-2.876652449688682) dist = 5.855266917542679e-07 normal = (0.8581505358745423,0.1418494641254577,0.4934129987196893) valid = 1
I found that GetLocalExitNormal missed rotating the normal from the 'next' volume (daughter-to-be, but not yet entered) to the 'current' volume (mother-to-be after the move). When I add it, the result now agrees with the Root result (for the case reported.)
Fix is included in tag geomnav-V10-04-04 and will be included in the next patch release.
