Ignore:
Timestamp:
Sep 10, 2008, 5:40:37 PM (17 years ago)
Author:
garnier
Message:

geant4.8.2 beta

Location:
trunk/source/geometry/solids/CSG
Files:
30 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/solids/CSG/History

    r831 r850  
    1 $Id: History,v 1.99.4.2 2008/04/23 09:55:26 gcosmo Exp $
     1$Id: History,v 1.105 2008/07/07 10:43:04 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1818     ----------------------------------------------------------
    1919
    20 Apr 23, 2008  G.Cosmo                    geom-csg-V09-00-03
    21 - Tag for release 9.1.p02.
     20Jul 07, 2008  V.Grichine                 geom-csg-V09-01-04
     21- G4Sphere: fixed bug in DistanceToOut(p, v, ...) for normal 'fSTheta'
     22  greater than 90*deg, and use of tangent giving negative value.
     23  Completes series of corrections introduced in previous tag.
     24
     25Jun 23, 2008  T.Nikitina                 geom-csg-V09-01-03
     26- G4Tubs: fix in DistanceToIn(p,v, ...) in case of point on surface with
     27  very small tangent direction; now returning kInfinity and no longer zero.
     28  It fixes rare observed cases of zero value returned by both DistanceToIn()
     29  and DistanceToOut(), causing stuck tracks with zero step length.
     30
     31Jun 20, 2008  V.Grichine
     32- G4Sphere: fixed calculation of roots in DistanceToOut(p,v,...) for
     33  theta-conical surfaces interserctions and for sTheta<=90 degree intersection.
     34  Addresses issue reported when running PET application with optical photons
     35  about mis-computation of distance on half-sphere constructs.
     36- Updated unit test for sphere.
    2237
    2338Apr 22, 2008  V.Grichine                 geom-csg-V09-01-02

  • trunk/source/geometry/solids/CSG/include/G4Box.hh

    r831 r850  
    2626//
    2727// $Id: G4Box.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Box.icc

    r831 r850  
    2626//
    2727// $Id: G4Box.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4CSGSolid.hh

    r831 r850  
    2626//
    2727// $Id: G4CSGSolid.hh,v 1.12 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// 

  • trunk/source/geometry/solids/CSG/include/G4Cons.hh

    r831 r850  
    2626//
    2727// $Id: G4Cons.hh,v 1.18 2007/05/18 07:38:00 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Cons.icc

    r831 r850  
    2626//
    2727// $Id: G4Cons.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Orb.hh

    r831 r850  
    2626//
    2727// $Id: G4Orb.hh,v 1.11 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Orb.icc

    r831 r850  
    2626//
    2727// $Id: G4Orb.icc,v 1.5 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Para.hh

    r831 r850  
    2626//
    2727// $Id: G4Para.hh,v 1.18 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Para.icc

    r831 r850  
    2626//
    2727// $Id: G4Para.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Sphere.hh

    r831 r850  
    2626//
    2727// $Id: G4Sphere.hh,v 1.20 2007/05/18 07:38:00 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Sphere.icc

    r831 r850  
    2626//
    2727// $Id: G4Sphere.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Torus.hh

    r831 r850  
    2626//
    2727// $Id: G4Torus.hh,v 1.27 2007/05/18 07:38:00 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Torus.icc

    r831 r850  
    2626//
    2727// $Id: G4Torus.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Trap.hh

    r831 r850  
    2626//
    2727// $Id: G4Trap.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Trap.icc

    r831 r850  
    2626//
    2727// $Id: G4Trap.icc,v 1.8 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Trd.hh

    r831 r850  
    2626//
    2727// $Id: G4Trd.hh,v 1.16 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Trd.icc

    r831 r850  
    2626//
    2727// $Id: G4Trd.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/include/G4Tubs.hh

    r831 r850  
    2626//
    2727// $Id: G4Tubs.hh,v 1.17 2007/05/18 07:38:00 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/include/G4Tubs.icc

    r831 r850  
    2626//
    2727// $Id: G4Tubs.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/src/G4Box.cc

    r831 r850  
    2626//
    2727// $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/src/G4CSGSolid.cc

    r831 r850  
    2626//
    2727// $Id: G4CSGSolid.cc,v 1.13 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/solids/CSG/src/G4Cons.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Cons.cc,v 1.54.4.1 2008/04/23 09:05:23 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4Cons.cc,v 1.56 2008/02/20 08:56:16 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/src/G4Orb.cc

    r831 r850  
    2525//
    2626// $Id: G4Orb.cc,v 1.24 2007/05/18 07:38:01 gcosmo Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// class G4Orb

  • trunk/source/geometry/solids/CSG/src/G4Para.cc

    r831 r850  
    2626//
    2727// $Id: G4Para.cc,v 1.39 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// class G4Para

  • trunk/source/geometry/solids/CSG/src/G4Sphere.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Sphere.cc,v 1.57.2.1 2008/04/23 09:05:23 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4Sphere.cc,v 1.68 2008/07/07 09:35:16 grichine Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// class G4Sphere
    … …  
    3434// History:
    3535//
     36// 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...)
    3637// 22.07.05 O.Link    : Added check for intersection with double cone
    3738// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
    … …  
    4748// 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections
    4849// 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections
    49 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
     50// 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
    5051// 17.09.96 V.Grichine: final modifications to commit
    5152// 28.03.94 P.Kent: old C++ code converted to tolerant geometry
    … …  
    608609      distETheta = std::fabs(pTheta-fSTheta-fDTheta);
    609610 
    610       nTs = G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
    611                         -std::cos(fSTheta)*std::sin(pPhi),
    612                          std::sin(fSTheta)               );
    613       nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
    614                          std::cos(fSTheta+fDTheta)*std::sin(pPhi),
    615                         -std::sin(fSTheta+fDTheta)               );   
     611      nTs = G4ThreeVector(-std::cos(fSTheta)*p.x()/rho, //  *std::cos(pPhi),
     612                          -std::cos(fSTheta)*p.y()/rho, //  *std::sin(pPhi),
     613                           std::sin(fSTheta)                   );
     614
     615      nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*p.x()/rho, // *std::cos(pPhi),
     616                           std::cos(fSTheta+fDTheta)*p.y()/rho, // *std::sin(pPhi),
     617                          -std::sin(fSTheta+fDTheta)                  );   
    616618    }
    617619    else if( !fRmin )
    618620    {
    619       if ( fSTheta )                distSTheta = 0.;
    620       if ( fSTheta + fDTheta < pi ) distETheta = 0.;
     621      if ( fSTheta ) 
     622      {             
     623        distSTheta = 0.;
     624        nTs = G4ThreeVector(0.,0.,-1.);
     625      }
     626      if ( fSTheta + fDTheta < pi ) // distETheta = 0.;
     627      {             
     628        distETheta = 0.;
     629        nTe = G4ThreeVector(0.,0.,1.);
     630      }
    621631    }   
    622632  }
    … …  
    14231433        d = std::sqrt(d2) ;
    14241434        s = -b - d ;    // First root
    1425 
    1426         if ( s < 0 )
    1427         {
    1428           s=-b+d;    // Second root
     1435        zi    = p.z() + s*v.z();
     1436
     1437        if ( s < 0 || zi*(fSTheta - halfpi) > 0 )
     1438        {
     1439          s = -b+d;    // Second root
    14291440        }
    14301441        if (s >= 0 && s < snxt)
    14311442        {
    1432           xi    = p.x() + s*v.x() ;
    1433           yi    = p.y() + s*v.y() ;
    1434           zi    = p.z() + s*v.z() ;
    1435           rhoi2 = xi*xi + yi*yi   ;
    1436           radi2 = rhoi2 + zi*zi   ;
     1443          xi    = p.x() + s*v.x();
     1444          yi    = p.y() + s*v.y();
     1445          zi    = p.z() + s*v.z();
     1446          rhoi2 = xi*xi + yi*yi;
     1447          radi2 = rhoi2 + zi*zi;
    14371448          if ( (radi2 <= tolORMax2)
    14381449            && (radi2 >= tolORMin2)
    … …  
    15011512      }
    15021513    } 
    1503     else if (pTheta > tolETheta)
    1504     { // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
    1505       // Inside (theta>etheta+tol) e theta cone
     1514    else if ( pTheta > tolETheta )
     1515    {
     1516      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
     1517      // Inside (theta > etheta+tol) e-theta cone
    15061518      // First root of etheta cone, second if first root `imaginary'
    15071519
    … …  
    15171529        d = std::sqrt(d2) ;
    15181530        s = -b - d ;    // First root
    1519         if (s < 0)
     1531        zi    = p.z() + s*v.z();
     1532
     1533        if (s < 0 || zi*(fSTheta + fDTheta - halfpi) > 0)
    15201534        {
    15211535          s = -b + d ;           // second root
    … …  
    19591973
    19601974  G4bool segTheta;                             // Theta flag and precals
    1961   G4double tanSTheta=0.,tanETheta, rhoSecTheta;
     1975  G4double tanSTheta=0.,tanETheta=0., rhoSecTheta;
    19621976  G4double tanSTheta2=0.,tanETheta2=0.;
    1963   G4double dist2STheta,dist2ETheta;
     1977  G4double dist2STheta, dist2ETheta, distTheta;
    19641978  G4double d2,s;
    19651979
    19661980  // General Precalcs
    19671981
    1968   rho2=p.x()*p.x()+p.y()*p.y();
    1969   rad2=rho2+p.z()*p.z();
     1982  rho2 = p.x()*p.x()+p.y()*p.y();
     1983  rad2 = rho2+p.z()*p.z();
    19701984  //  G4double rad=std::sqrt(rad2);
    19711985
    1972   pTheta=std::atan2(std::sqrt(rho2),p.z());
    1973 
    1974   pDotV2d=p.x()*v.x()+p.y()*v.y();
    1975   pDotV3d=pDotV2d+p.z()*v.z();
     1986  pTheta = std::atan2(std::sqrt(rho2),p.z());
     1987
     1988  pDotV2d = p.x()*v.x()+p.y()*v.y();
     1989  pDotV3d = pDotV2d+p.z()*v.z();
    19761990
    19771991  // Set phi divided flag and precalcs
    19781992
    1979   if(fDPhi<twopi)
     1993  if( fDPhi < twopi )
    19801994  {
    19811995    segPhi=true;
    … …  
    19962010  // Theta precalcs
    19972011   
    1998   if (fDTheta < pi)
    1999   {
    2000     segTheta=true;
    2001     tolSTheta=fSTheta-kAngTolerance*0.5;
    2002     tolETheta=fSTheta+fDTheta+kAngTolerance*0.5;
    2003   }
    2004   else
    2005   {
    2006     segTheta=false;
    2007   }
     2012  if ( fDTheta < pi )
     2013  {
     2014    segTheta  = true;
     2015    tolSTheta = fSTheta - kAngTolerance*0.5;
     2016    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5;
     2017  }
     2018  else segTheta = false;
     2019
    20082020   
    20092021  // Radial Intersections from G4Sphere::DistanceToIn
    … …  
    20242036  //
    20252037  // const G4double  fractionTolerance = 1.0e-12;
    2026   const G4double  flexRadMaxTolerance = // kRadTolerance;
     2038
     2039  const G4double  flexRadMaxTolerance =            // kRadTolerance;
    20272040    std::max(kRadTolerance, fEpsilon * fRmax);
    20282041
    20292042  const G4double  Rmax_plus = fRmax + flexRadMaxTolerance*0.5;
     2043
    20302044  const G4double  flexRadMinTolerance = std::max(kRadTolerance,
    20312045                     fEpsilon * fRmin);
     2046
    20322047  const G4double  Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;
    20332048
    … …  
    20632078      else
    20642079      {
    2065         snxt=-pDotV3d+std::sqrt(d2);    // second root since inside Rmax
    2066         side = kRMax ;
     2080        snxt = -pDotV3d+std::sqrt(d2);    // second root since inside Rmax
     2081        side =  kRMax ;
    20672082      }
    20682083    }
    … …  
    20772092      d2 = pDotV3d*pDotV3d - c;
    20782093
    2079       if (c >- flexRadMinTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
     2094      if ( c >- flexRadMinTolerance*fRmin ) // 2.0 * (0.5*kRadTolerance) * fRmin
    20802095      {
    20812096        if( c < flexRadMinTolerance*fRmin &&
    20822097            d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
    20832098        {
    2084           if(calcNorm)
    2085           {
    2086             *validNorm = false ;   // Rmin surface is concave
    2087           }
    2088           return snxt = 0 ;
     2099          if(calcNorm)  *validNorm = false ;   // Rmin surface is concave         
     2100                       return snxt = 0 ;
    20892101        }
    20902102        else
    20912103        { 
    2092           if (d2 >= 0)
    2093           {
    2094             s = -pDotV3d-std::sqrt(d2) ;
    2095             if (s>=0)     // Always intersect Rmin first
     2104          if ( d2 >= 0. )
     2105          {
     2106            s = -pDotV3d-std::sqrt(d2);
     2107
     2108            if ( s >= 0. )     // Always intersect Rmin first
    20962109            {
    20972110              snxt = s ;
    … …  
    21292142    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
    21302143    //
     2144 
     2145    /* ////////////////////////////////////////////////////////
     2146
    21312147    tanSTheta=std::tan(fSTheta);
    21322148    tanSTheta2=tanSTheta*tanSTheta;
    … …  
    22882304              s = kInfinity ;  // wrong cone
    22892305            }
    2290             if (s < stheta)
    2291             {
    2292               stheta = s ;
    2293               sidetheta = kETheta ;
    2294             }
     2306          }
     2307          if (s < stheta)
     2308          {
     2309            stheta = s ;
     2310            sidetheta = kETheta ;
    22952311          }
    22962312        }
    22972313      }
    22982314    } 
    2299   }
     2315    */  ////////////////////////////////////////////////////////////
     2316
     2317    if(fSTheta) // intersection with first cons
     2318    {
     2319
     2320      tanSTheta = std::tan(fSTheta);
     2321
     2322      if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
     2323      {
     2324        if( v.z() > 0. )
     2325        {
     2326          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
     2327          {
     2328            if(calcNorm)
     2329            {
     2330              *validNorm = true;
     2331              *n = G4ThreeVector(0.,0.,1.);
     2332            }
     2333            return snxt = 0 ;
     2334          } 
     2335          // s = -p.z()/v.z();
     2336          stheta    = -p.z()/v.z();
     2337          sidetheta = kSTheta;
     2338        }
     2339      }
     2340      else // kons is not plane
     2341      {
     2342        tanSTheta2  = tanSTheta*tanSTheta;
     2343        t1          = 1-v.z()*v.z()*(1+tanSTheta2);
     2344        t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
     2345        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;      // t3
     2346
     2347        // distTheta = std::sqrt(std::fabs(dist2STheta/(1+tanSTheta2)));
     2348        distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
     2349
     2350        if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
     2351        {
     2352          if( v.z() > 0. )
     2353          {
     2354            if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
     2355            {
     2356              if( fSTheta < halfpi && p.z() > 0. )
     2357              {
     2358                if( calcNorm ) *validNorm = false;
     2359                              return snxt = 0.;
     2360              }
     2361              else if( fSTheta > halfpi && p.z() <= 0)
     2362              {
     2363                if( calcNorm )
     2364                {
     2365                  *validNorm = true;
     2366                  if (rho2)
     2367                  {
     2368                    rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
     2369                   
     2370                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
     2371                                        p.y()/rhoSecTheta,
     2372                                        std::sin(fSTheta)  );
     2373                  }
     2374                  else *n = G4ThreeVector(0.,0.,1.);
     2375                }
     2376                return snxt = 0.;               
     2377              }
     2378            }
     2379            // s = -0.5*dist2STheta/t2;
     2380
     2381            stheta    = -0.5*dist2STheta/t2;
     2382            sidetheta = kSTheta;
     2383          } 
     2384        }
     2385        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
     2386        {
     2387          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
     2388          {
     2389            if( fSTheta > halfpi && t2 >= 0. ) // leave
     2390            {
     2391              if( calcNorm )
     2392              {
     2393                *validNorm = true;
     2394                if (rho2)
     2395                {
     2396                    rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
     2397                   
     2398                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
     2399                                        p.y()/rhoSecTheta,
     2400                                        std::sin(fSTheta)  );
     2401                }
     2402                else *n = G4ThreeVector(0.,0.,1.);
     2403              }                           
     2404              return snxt = 0.;
     2405            }
     2406            else if( fSTheta < halfpi  && t2 < 0. && p.z() >=0. ) // leave
     2407            {
     2408                if( calcNorm )   *validNorm = false;                                                 
     2409                return snxt = 0.;
     2410            }                               
     2411          }
     2412          b  = t2/t1;
     2413          c  = dist2STheta/t1;
     2414          d2 = b*b - c ;
     2415
     2416          if ( d2 >= 0. )
     2417          {
     2418            d = std::sqrt(d2);
     2419
     2420            if( fSTheta > halfpi )
     2421            {
     2422              s = -b - d;         // First root
     2423
     2424              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
     2425                              s < 0.  ||
     2426                  ( s > 0. && p.z() + s*v.z() > 0.)                   )
     2427              {
     2428                s = -b + d ; // 2nd root
     2429              }
     2430              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 
     2431              {
     2432                stheta    = s;
     2433                sidetheta = kSTheta;
     2434              }
     2435            }
     2436            else // sTheta < pi/2, concave surface, no normal
     2437            {
     2438              s = -b - d;         // First root
     2439
     2440              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
     2441                              s < 0.                                   ||
     2442                  ( s > 0. && p.z() + s*v.z() < 0.)                      )
     2443              {
     2444                s = -b + d ; // 2nd root
     2445              }
     2446              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() >= 0.) 
     2447              {
     2448                stheta    = s;
     2449                sidetheta = kSTheta;
     2450              }           
     2451            }
     2452          }
     2453        }
     2454      }
     2455    }
     2456    if (fSTheta + fDTheta < pi) // intersection with second cons
     2457    {
     2458
     2459      tanETheta = std::tan(fSTheta+fDTheta);
     2460
     2461      if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
     2462      {
     2463        if( v.z() < 0. )
     2464        {
     2465          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
     2466          {
     2467            if(calcNorm)
     2468            {
     2469              *validNorm = true;
     2470              *n = G4ThreeVector(0.,0.,-1.);
     2471            }
     2472            return snxt = 0 ;
     2473          } 
     2474          s = -p.z()/v.z();
     2475
     2476          if( s < stheta)
     2477          {
     2478            stheta    = s;
     2479            sidetheta = kETheta;
     2480          }
     2481        }
     2482      }
     2483      else // kons is not plane
     2484      {
     2485        tanETheta2  = tanETheta*tanETheta;
     2486        t1          = 1-v.z()*v.z()*(1+tanETheta2);
     2487        t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
     2488        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;      // t3
     2489
     2490        // distTheta = std::sqrt(std::fabs(dist2ETheta/(1+tanETheta2)));
     2491        distTheta = std::sqrt(rho2)-p.z()*tanETheta;
     2492
     2493        if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
     2494        {
     2495          if( v.z() < 0. )
     2496          {
     2497            if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
     2498            {
     2499              if( fSTheta+fDTheta > halfpi && p.z() < 0. )
     2500              {
     2501                if( calcNorm ) *validNorm = false;
     2502                              return snxt = 0.;
     2503              }
     2504              else if( fSTheta+fDTheta < halfpi && p.z() >= 0)
     2505              {
     2506                if( calcNorm )
     2507                {
     2508                  *validNorm = true;
     2509                  if (rho2)
     2510                  {
     2511                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
     2512                   
     2513                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
     2514                                        p.y()/rhoSecTheta,
     2515                                        -std::sin(fSTheta+fDTheta)  );
     2516                  }
     2517                  else *n = G4ThreeVector(0.,0.,-1.);
     2518                }
     2519                return snxt = 0.;               
     2520              }
     2521            }
     2522            s = -0.5*dist2ETheta/t2;
     2523
     2524            if( s < stheta)
     2525            {
     2526              stheta    = s;
     2527              sidetheta = kETheta;
     2528            }
     2529          } 
     2530        }
     2531        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
     2532        {
     2533          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
     2534          {
     2535            if( fSTheta+fDTheta < halfpi && t2 >= 0. ) // leave
     2536            {
     2537              if( calcNorm )
     2538              {
     2539                *validNorm = true;
     2540                if (rho2)
     2541                {
     2542                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
     2543                   
     2544                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
     2545                                        p.y()/rhoSecTheta,
     2546                                        -std::sin(fSTheta+fDTheta)  );
     2547                }
     2548                else *n = G4ThreeVector(0.,0.,-1.);
     2549              }                           
     2550              return snxt = 0.;
     2551            }
     2552            else if( fSTheta+fDTheta > halfpi  && t2 < 0. && p.z() <=0. ) // leave
     2553            {
     2554                if( calcNorm )   *validNorm = false;                                                 
     2555                return snxt = 0.;
     2556            }                               
     2557          }
     2558          b  = t2/t1;
     2559          c  = dist2ETheta/t1;
     2560          d2 = b*b - c ;
     2561
     2562          if ( d2 >= 0. )
     2563          {
     2564            d = std::sqrt(d2);
     2565
     2566            if( fSTheta+fDTheta < halfpi )
     2567            {
     2568              s = -b - d;         // First root
     2569
     2570              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
     2571                              s < 0. )
     2572              {
     2573                s = -b + d ; // 2nd root
     2574              }
     2575              if( s >  flexRadMaxTolerance*0.5 ) 
     2576              {
     2577                if( s < stheta )
     2578                {
     2579                  stheta    = s;
     2580                  sidetheta = kETheta;
     2581                }
     2582              }
     2583            }
     2584            else // sTheta+fDTheta > pi/2, concave surface, no normal
     2585            {
     2586              s = -b - d;         // First root
     2587
     2588              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
     2589                              s < 0.                                   ||
     2590                  ( s > 0. && p.z() + s*v.z() > 0.)                      )
     2591              {
     2592                s = -b + d ; // 2nd root
     2593              }
     2594              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 
     2595              {
     2596                if( s < stheta )
     2597                {
     2598                  stheta    = s;
     2599                  sidetheta = kETheta;
     2600                }
     2601              }           
     2602            }
     2603          }
     2604        }
     2605      }
     2606    }
     2607
     2608  } // end theta intersections
    23002609
    23012610  // Phi Intersection
    … …  
    25782887        *validNorm=true;
    25792888        break;
     2889
    25802890      case kRMin:
    25812891        *validNorm=false;  // Rmin is concave
    25822892        break;
     2893
    25832894      case kSPhi:
    2584         if (fDPhi<=pi)     // Normal to Phi-
     2895        if ( fDPhi <= pi )     // Normal to Phi-
    25852896        {
    25862897          *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
    … …  
    25892900        else *validNorm=false;
    25902901        break ;
     2902
    25912903      case kEPhi:
    2592         if (fDPhi<=pi)      // Normal to Phi+
     2904        if ( fDPhi <= pi )      // Normal to Phi+
    25932905        {
    25942906          *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
    … …  
    25972909        else *validNorm=false;
    25982910        break;
     2911
    25992912      case kSTheta:
    2600         if( fSTheta == pi*0.5 )
    2601         {
    2602           *n=G4ThreeVector(0,0,1);
     2913        if( fSTheta == halfpi )
     2914        {
     2915          *n=G4ThreeVector(0.,0.,1.);
    26032916          *validNorm=true;
    26042917        }
    2605         else if ( fSTheta > pi )
     2918        else if ( fSTheta > halfpi )
     2919        {
     2920          xi = p.x() + snxt*v.x();
     2921          yi = p.y() + snxt*v.y();
     2922          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2));
     2923          *n = G4ThreeVector( xi/rhoSecTheta,   // N-
     2924                              yi/rhoSecTheta,
     2925                              -tanSTheta/std::sqrt(1+tanSTheta2));
     2926          *validNorm=true;
     2927        }
     2928        else *validNorm=false;  // Concave STheta cone
     2929        break;
     2930
     2931      case kETheta:
     2932        if( ( fSTheta + fDTheta ) == halfpi )
     2933        {
     2934          *n         = G4ThreeVector(0.,0.,-1.);
     2935          *validNorm = true;
     2936        }
     2937        else if ( ( fSTheta + fDTheta ) < halfpi)
    26062938        {
    26072939          xi=p.x()+snxt*v.x();
    26082940          yi=p.y()+snxt*v.y();
    2609           rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2)) ;
    2610           *n = G4ThreeVector(-xi/rhoSecTheta,   // N-
    2611                              -yi/rhoSecTheta,
    2612                              tanSTheta/std::sqrt(1+tanSTheta2)) ;
    2613           *validNorm=true;
    2614         }
    2615         else *validNorm=false;  // Concave STheta cone
    2616         break;
    2617       case kETheta:
    2618         if( ( fSTheta + fDTheta ) == pi*0.5 )
    2619         {
    2620           *n         = G4ThreeVector(0,0,-1);
    2621           *validNorm = true ;
    2622         }
    2623         else if ( ( fSTheta + fDTheta ) < pi )
    2624         {
    2625           xi=p.x()+snxt*v.x();
    2626           yi=p.y()+snxt*v.y();
    2627           rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2)) ;
     2941          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2));
    26282942          *n = G4ThreeVector( xi/rhoSecTheta,   // N+
    26292943                              yi/rhoSecTheta,
    2630                               -tanSTheta/std::sqrt(1+tanSTheta2) ) ;
     2944                              -tanETheta/std::sqrt(1+tanETheta2) );
    26312945          *validNorm=true;
    26322946        }
    26332947        else *validNorm=false;   // Concave ETheta cone
    26342948        break;
     2949
    26352950      default:
    26362951        G4cout.precision(16);

  • trunk/source/geometry/solids/CSG/src/G4Torus.cc

    r831 r850  
    2626//
    2727// $Id: G4Torus.cc,v 1.63 2007/10/02 09:34:17 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/src/G4Trap.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Trap.cc,v 1.42.4.1 2008/04/23 09:55:26 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4Trap.cc,v 1.45 2008/04/23 09:49:57 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// class G4Trap

  • trunk/source/geometry/solids/CSG/src/G4Trd.cc

    r831 r850  
    2626//
    2727// $Id: G4Trd.cc,v 1.34 2006/10/19 15:33:38 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/CSG/src/G4Tubs.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Tubs.cc,v 1.66 2007/11/23 09:07:43 tnikitin Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4Tubs.cc,v 1.68 2008/06/23 13:37:39 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    983983          iden   = std::sqrt(t3) ;
    984984          cosPsi = inum/iden ;
    985           if (cosPsi >= cosHDPhiIT)  { return 0.0; }
     985          if (cosPsi >= cosHDPhiIT)
     986          {
     987            // In the old version, the small negative tangent for the point
     988            // on surface was not taken in account, and returning 0.0 ...
     989            // New version: check the tangent for the point on surface and
     990            // if no intersection, return kInfinity, if intersection instead
     991            // return s.
     992            //
     993            c = t3-fRMax*fRMax;
     994            if ( c<=0.0 )
     995            {
     996              return 0.0;
     997            }
     998            else
     999            {
     1000              c = c/t1 ;
     1001              d = b*b-c;
     1002              if ( d>=0.0 )
     1003              {
     1004                snxt = c/(-b+std::sqrt(d)); // using safe solution
     1005                                            // for quadratic equation
     1006                if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
     1007                return snxt ;
     1008              }     
     1009              else
     1010              {
     1011                return kInfinity;
     1012              }
     1013            }
     1014          }
    9861015        }
    9871016        else
    988         {
    989           return 0.0 ;
    990         }
    991       }
    992     }     
     1017        {   
     1018          // In the old version, the small negative tangent for the point
     1019          // on surface was not taken in account, and returning 0.0 ...
     1020          // New version: check the tangent for the point on surface and
     1021          // if no intersection, return kInfinity, if intersection instead
     1022          // return s.
     1023          //
     1024          c = t3 - fRMax*fRMax;
     1025          if ( c<=0.0 )
     1026          {
     1027            return 0.0;
     1028          }
     1029          else
     1030          {
     1031            c = c/t1 ;
     1032            d = b*b-c;
     1033            if ( d>=0.0 )
     1034            {
     1035              snxt= c/(-b+std::sqrt(d)); // using safe solution
     1036                                         // for quadratic equation
     1037              if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
     1038              return snxt ;
     1039            }     
     1040            else
     1041            {
     1042              return kInfinity;
     1043            }
     1044          }
     1045        } // end if   (seg)
     1046      }   // end if   (t3>tolIRMin2)
     1047    }     // end if   (Inside Outer Radius)
    9931048    if ( fRMin )    // Try inner cylinder intersection
    9941049    {
    … …  
    14591514              //
    14601515              if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
    1461                 { sidephi = kSPhi;
     1516                { sidephi = kSPhi;
    14621517                if (((fSPhi-0.5*kAngTolerance)<=vphi)
    14631518                   &&((ePhi+0.5*kAngTolerance)>=vphi))
Note: See TracChangeset for help on using the changeset viewer.