The point is (0,0,0), the solid is the following: <Polycone name="ESPM" startPhi="-10.25*deg" deltaPhi="21*deg" > <ZSection z="-1.52*m" rMin="1.238*m" rMax="1.55501*m" /> <ZSection z="-80.4*cm" rMin="1.238*m" rMax="1.55501*m" /> <ZSection z="-80.4*cm" rMin="1.238*m" rMax="1.53805*m" /> <ZSection z="-51.5345*cm" rMin="1.238*m" rMax="1.53805*m" /> <ZSection z="-51.5345*cm" rMin="1.238*m" rMax="1.52326*m" /> <ZSection z="-17.7*cm" rMin="1.238*m" rMax="1.52326*m" /> <ZSection z="-17.7*cm" rMin="1.238*m" rMax="1.50624*m" /> <ZSection z="14.9561*cm" rMin="1.238*m" rMax="1.50624*m" /> <ZSection z="14.9561*cm" rMin="1.238*m" rMax="1.48852*m" /> <ZSection z="57.5*cm" rMin="1.238*m" rMax="1.48852*m" /> <ZSection z="57.5*cm" rMin="1.238*m" rMax="1.47128*m" /> <ZSection z="98.2812*cm" rMin="1.238*m" rMax="1.47128*m" /> <ZSection z="98.2812*cm" rMin="1.238*m" rMax="1.45522*m" /> <ZSection z="1.1667*m" rMin="1.238*m" rMax="1.45522*m" /> <ZSection z="1.524*m" rMin="1.45522*m" rMax="1.45522*m" /> </Polycone> It is placed at <Translation x="0*fm" y="0*fm" z="1.52*m" /> with no rotation I have put some print outs in G4PolyconeSide::DistanceAway to understand what happens. Here it is the modified method G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p, G4bool opposite, G4double &distOutside2, G4double *edgeRZnorm ) { int myverbose = 0; char* g4evtv = getenv("G4EVENTVERBOSE"); int ig4evtverb = atoi(g4evtv); if( G4EventManager::GetEventManager() ) { if( G4EventManager::GetEventManager()->GetConstCurrentEvent() ) { if( G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()> ig4evtverb ) myverbose = 2; } } // // Convert our point to r and z // G4double rx = p.perp(), zx = p.z(); // // Change sign of r if opposite says we should // if (opposite) rx = -rx; // // Calculate return value // G4double deltaR = rx - r[0], deltaZ = zx - z[0]; G4double answer = deltaR*rNorm + deltaZ*zNorm; // // Are we off the surface in r,z space? // G4double s = deltaR*rS + deltaZ*zS; if( myverbose >= 1 ) G4cout << p << " " << r[0] << " " << z[0] << " DistanceAway s " << s << " " << deltaR << " " << rS << " " << deltaZ << " " << zS << G4endl; if (s < 0){ distOutside2 = s*s; if( myverbose >= 1 ) G4cout << " s < 0 " << s << " " << length << " " << distOutside2 << G4endl; if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0]; }else if (s > length) { distOutside2 = sqr( s-length ); if( myverbose >= 1 )G4cout << " s > length " << s << " " << length << " " << distOutside2 << G4endl; if (edgeRZnorm){ G4double deltaR = rx - r[1], deltaZ = zx - z[1]; *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1]; } } else{ distOutside2 = 0; if( myverbose >= 1 )G4cout << " s < length " << s << " " << length << " " << distOutside2 << G4endl; if (edgeRZnorm) *edgeRZnorm = answer; } if (phiIsOpen){ // // Finally, check phi // G4double phi = p.phi(); while( phi < startPhi ) phi += 2*M_PI; if (phi > startPhi+deltaPhi) { // // Oops. Are we closer to the start phi or end phi? // G4double d1 = phi-startPhi-deltaPhi; while( phi > startPhi ) phi -= 2*M_PI; G4double d2 = startPhi-phi; if (d2 < d1) d1 = d2; // // Add result to our distance // G4double dist = d1*rx; distOutside2 += dist*dist; if (edgeRZnorm) *edgeRZnorm = fabs(dist); } } return answer; } ---------- And this is the relevant output: (-9.30701e-14,-9.30701e-14,-1520) 1455.22 1524 DistanceAway s 3357 -1455.22 -0.519481 -3044 -0.854482 s > length 3357 418.148 8.63686e+06 (-9.30701e-14,-9.30701e-14,-1520) 1238 1166.7 DistanceAway s 2686.7 -1238 0 -2686.7 -1 s < length 2686.7 2686.7 0 It seems that while checking the second face (the one at rmin = 1238 and between z= -1520, z = 1166.7, it finds the correct distance in R and Z but when it multiplies them by the vector normal to the surface, that is (0,1), the resulting distance is equal to the face length, what makes distOutside2 = 0 and then it thinks it is on the surface
There could be a problem, but the reporter's description does not correctly indicate any. His attempts to diagnosis the code, however, are appreciated. The variable "distOutside2" is the distance away from the surface perpendicular to the surface normal. The quantity of greater interest is the return value of the function itself (stored in the variable "answer"), which is the distance of the point from the surface along the normal of the surface. I have attempted to reproduce the given geometry and point, and in my code the point is correctly returned as "outside." It could be that the issue is related to floating point roundoff, which is why it might be difficult to reproduce. I will contact the error reporter directly and request additional information.
I believe I have found this problem and have fixed it in CVS (revision 1.7 of geometry/solids/specific/ src/G4PolyconeSide.cc).
F.V.Ignatov@inp.nsk.su writes: With the new geant4.7.0 G4Polycone::Inside() still returns kSurface when an input point is outside and on the end of "phi angle". Example of incorrect behaviour: #include "G4Polycone.hh" #include "G4ThreeVector.hh" #include <iostream> int main() { double z[2]={0,100}; double ri[2]={510.,510.}; double ro[2]={519.,519.}; G4Polycone *solidStore = new G4Polycone("TEST",(-360./33*0.5)*deg, (360./33*deg),2,z,ri,ro); G4ThreeVector point(504.116945822705,-48.1373317980017,50); std::cout<<*solidStore<<std::endl; std::cout<<point<<std::endl; std::cout<<(kSurface==solidStore->Inside(point))<<" " <<(kInside==solidStore->Inside(point))<<std::endl; exit(1); } make output: ----------------------------------------------------------- *** Dump for solid - TEST *** =================================================== Solid type: G4Polycone Parameters: starting phi angle : 354.545 degrees ending phi angle : 365.455 degrees number of Z planes: 2 Z values: Z plane 0: 0 Z plane 1: 100 Tangent distances to inner surface (Rmin): Z plane 0: 510 Z plane 1: 510 Tangent distances to outer surface (Rmax): Z plane 0: 519 Z plane 1: 519 number of RZ points: 4 RZ values (corners): 510, 100 510, 0 519, 0 519, 100 ----------------------------------------------------------- (504.117,-48.1373,50) 1 0 I think the last correction before geant4.7.0 in G4PolyconeSide::DistanceAway() doesn't resolve problem: if (edgeRZnorm) *edgeRZnorm = std::max(*edgeRZnorm,std::fabs(dist)); if "*edgeRZnorm" is negative then the result will be zero, even if point is not at G4PolyconeSide, so it must be like this: if (edgeRZnorm) *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist)); Sincerely yours, Fedor Ignatov
I want to thank Pedro for supplying the example code. This is always very helpful. Unfortunately, when I compile this example on my system, I cannot reproduce the error. This makes the problem a little more difficult to diagnose. However, after studying the issue, I think I understand the problem. First, let me talk about your suggestion. The proposed code change: if (edgeRZnorm) *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist)); is not correct, as edgeRZnorm is a signed parametered, and the sign has important significance. If this code change was made, it will break many things. Now, back to the issue. The way G4Polycone works is that the solid surface that is closest to the point is used to decide whether the point is inside or outside the surface. To decide which surface is closest, a "distance" value is calculated. This distance is meant to be the closest point the solid surface approaches the point, and thus needs to account for the edges of each surface. You will note that the point in your example is equally close to two surfaces, thus, depending on floating point roundoff, either surface might be chosen to decide if the point is inside or not. The problem is that the point is actually just "inside" one of these surfaces, since inside/outside is logically defined in reference to the extrapolation of each surface in all directions. I apologize if this explanation does not seem sufficient. Needless to say, it is a serious design flaw in G4VCSGfaceted. I do not know right now how to fix it, but I continue to work on this. Again, thank you for reporting this bug.
The problem is fixed and submitted in tag "geom-specific-V08-02-01". It will be available with the next release. Bug was found in G4PolyconeSide::DistanceAway in case of Open Polycone in phi: if (edgeRZnorm) *edgeRZnorm = std::max(*edgeRZnorm,std::fabs(dist)); In case of OpenPolycone the 'edgeRZnorm' could be wrong calculated. Example: Point situated outside of the Polycone in R and Z directions, but close to the Polycone in phi direction. When âInsideâ function is called, the loop on all faces is performed. If 'edgeRZnorm'is is negative for one of Polycone faces, according the algorithm 'edgeRZnormâ will return |dist in phi|. If this distance is very small, the 'Inside' function will return 'kSurface' when the point is outside. To fix this bug the proposal of F.V.Ignatov is correct : if (edgeRZnorm) *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs (dist)); It can be used, because we can arrive to this line in the algorithm only if the point is outside of the Polycone in phi. This point can be only outside or on the surface of the Polycone. Since positive 'edgeRZnorm' means outside, positive values of 'edgeRZnorm' can be used. The calculated value of 'edgeRZnorm' is used only in 'Inside' function and gives a good answer if point is outside or on the surface. If the point is not on the surface, the distance to phi-edge is calculated using G4PolyPhiFace.cc G4BREPSolidPcone and G4Polyhedra were checked, they donât have this mistake in their algorithm.