Problem 1025

The bug comes from G4TriangularFacet::GetPointOnFace() method. Here is the buggy implementation (from G4TriangularFacet.cc):

---------------------------------
G4ThreeVector G4TriangularFacet::GetPointOnFace() const   
{
  G4double lambda0 = CLHEP::RandFlat::shoot(0.,1.);
  G4double lambda1 = CLHEP::RandFlat::shoot(0.,lambda0);

  return (P0 + lambda0*E[0] + lambda1*E[1]);
}
---------------------------------

This class assumes a triangle defined by three points (P0, P1, P2), or a point P0 and two vectors E[0]=(P1-P0) and E[1]=(P2-P0), using the notation from the implementation above. The function should return a point on the triangle or inside it. The parameters lambda0,lambda1 are both taken randomly. Take for instance some trivial cases for lambda0 and lambda1 parameters:

    * lambda0=0, lambda1=0 gives vertex P0
    * lambda0=1, lambda1=0 gives vertex P1 = P0+E[0]
    * P2=P0+E[1] is never obtained, because lambda0=0 and lambda1=1 cannot be generated by the code above
    * lambda0=1, lambda1=1 returns a point outside the triangle: P0 + E[0] + E[1] (wrong!) 


With points returned outside the triangle, geometry checks like G4PVPlacement::CheckOverlaps() may return many false positives for solids based on facets (like G4ExtrudedSolids and similars).

=================================

I have two proposals to fix this bug:

    * Option 1: 

G4ThreeVector G4TriangularFacet::GetPointOnFace() const
{
  G4double lambda0 = CLHEP::RandFlat::shoot(0.,1.);
  G4double lambda1 = CLHEP::RandFlat::shoot(0.,1-lambda0);

  return (P0 + lambda0*E[0] + lambda1*E[1]);
}

    * Option 2: 

G4ThreeVector G4TriangularFacet::GetPointOnFace() const
{
  G4double alpha = CLHEP::RandFlat::shoot(0.,1.);
  G4double beta  = CLHEP::RandFlat::shoot(0.,1.);
  G4double lambda1 = alpha*beta;
  G4double lambda0 = alpha - lambda1;

  return (P0 + lambda0*E[0] + lambda1*E[1]);
}


=======================================================

Derivation of the bug fix
-------------------------

In order to fix the bug, one has to somehow introduce the line P1,P2 as a constraint. The vector equation for the (infinite) line going through points P1,P2, in terms of vectors P0,E[0],E[1], is this:

   (P0 + E[0]) + beta*( E[1]-E[0] )

since P0+E[0] is P1 (point P1 belongs to line) and E[1]-E[0] is the vector defining the line orientation in space. Parameter beta is any real number for the infinite line, but as shown above, beta has to be in the range [0,1] to span the points between P1 and P2. Rewritting the expression above one gets:

   P0 + (1-beta)*E[0] + beta*E[1]

representing the line from P1 to P2 with respect to P0.

To span the triangle one can just introduce another parameter alpha, also in the range [0,1], which brings one from P0 (alpha=0) to the line (alpha=1):

   P0 + alpha*( (1-beta)*E[0] + beta*E[1] )

or, as suggested in the code:

   P0 + lambda0*E[0] + lambda1*e[1]

where:

   lambda0 = alpha*(1-beta)
   lambda1 = alpha*beta,

or

lambda0 = alpha - lambda1

(this is the triangle restriction). These expressions already suggests the second algorithm option as a bug fix.

But I actually found option (1) first than option (2), and this is how I found it: the parameters lambda0,lambda1 are randoms, therefore all one needs to know is their validity interval. Replacing alpha above is tricky, because lambda0,lambda1 cannot be negative.

I found easier to look at the "inside triangle" restriction as 0<alpha<1, for determining lambda0,lambda1 ranges:

   0 < lambda0 < 1-beta
   0 < lambda1 < beta

or if lambda1=1, lambda0 can only be zero. So here is the first algorithm: draw lambda1 (beta) in the range [0,1], and then draw lambda0 in the range [0,1-lambda1]. That's the algorithm option 1.


Thanks,
Guilherme