Problem 1073 - Special case not managed in Polycone side
Summary: Special case not managed in Polycone side
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: geometry/solids (show other problems)
Version: 9.0
Hardware: All All
: P5 trivial
Assignee: tatiana.nikitina
URL:
Depends on:
Blocks:
 
Reported: 2009-07-03 23:06 CEST by Pierre
Modified: 2009-07-24 17:32 CEST (History)
1 user (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description Pierre 2009-07-03 23:06:49 CEST
Hello,
this bug occurs in a very particular situation when the Polycone using the plane (Z=0) as an extremity. In that case, the normal at the center of the extrimity is given by the function :

G4ThreeVector G4PolyconeSide::Normal( const G4ThreeVector &p,
                                            G4double *bestDistance )
{
  if (p == G4ThreeVector(0.,0.,0.))  { return p; }

  G4ThreeVector dFrom;
  G4double dOut2;

  dFrom = DistanceAway( p, false, dOut2 );

  *bestDistance = std::sqrt( dFrom*dFrom + dOut2 );

  G4double rad = p.perp();
  return G4ThreeVector( rNorm*p.x()/rad, rNorm*p.y()/rad, zNorm );
}

with p=(0,0,0) so rad=(0,0,0) and the Vector return is (nan,nan,0). 

I changed it to 

  if (rad != 0)
    return G4ThreeVector( rNorm*p.x()/rad, rNorm*p.y()/rad, zNorm );
  else
    return G4ThreeVector( rNorm*p.x()*1000000., rNorm*p.y()*1000000., zNorm );

which gives acceptable results, but you may come up with a better solution.

Regards,
Pierre
Comment 1 tnikitin 2009-07-21 16:18:21 CEST
Dear Pierre,

 Thank you for reporting this problem.

 For p(0,0,z) rad=p.perp() is equal to zero.
 In calculation of Normal division by rad is causing nan:
  return G4ThreeVector( rNorm*p.x()/rad, rNorm*p.y()/rad, zNorm ).

 But in this case p.x() and p.y() are zeros and Normal can be given by :
  return G4ThreeVector( 0,0, zNorm ).unit().

 Fix will be provided in the next Geant4 release.

 Best Regards,
 
 Tatiana Nikitina