Problem 1089

Greetings,
In G4Ellipsoid, a recent modification was made to correct, I believe, the bug #1050. I have some reserve about the line 522

if (p.dot(v)<0.) { distMin = distR; }

it does not make sense to compare p and v but rather to compare v and the normal to the surface in order to know if we are going toward it or not...

I would then suggest 

  G4double A,B,C;

  A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
  C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
  B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
           + p.y()*v.y()/(ySemiAxis*ySemiAxis)
           + p.z()*v.z()/(zSemiAxis*zSemiAxis) );

  C= B*B - 4.0*A*C;
  if (C > 0.0)
  {    
    G4double distR= (-B - std::sqrt(C)) / (2.0*A);
    G4double intZ = p.z()+distR*v.z();
    if ( (distR >- halfRadTolerance)
      && (intZ >= zBottomCut-halfRadTolerance)
      && (intZ <= zTopCut+halfRadTolerance) )
    { 
      distMin = distR;
    }
    else
    {
      distR= (-B + std::sqrt(C)) / (2.0*A);
      intZ = p.z()+distR*v.z();
      if ( (distR >- halfRadTolerance)
        && (intZ >= zBottomCut-halfRadTolerance)
        && (intZ <= zTopCut+halfRadTolerance) )
      {
        G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
                           p.y()/(ySemiAxis*ySemiAxis),
                           p.z()/(zSemiAxis*zSemiAxis));
        if (norm.dot(v)<0.) {  { distMin = distR; }
      }
    }
    if ( distMin>dRmax ) // Avoid rounding errors due to precision issues on
    {                    // 64 bits systems. Split long distances and recompute
      G4double fTerm = distMin-std::fmod(distMin,dRmax);
      distMin = fTerm + DistanceToIn(p+fTerm*v,v);
    }
  }

Best regards,
Pierre