Problem 598

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

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.