Ignore:
Timestamp:
Feb 16, 2009, 10:14:30 AM (15 years ago)
Author:
garnier
Message:

en test de gl2ps. Problemes de libraries

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r850 r921  
    2525//
    2626//
    27 // $Id: G4Tubs.cc,v 1.68 2008/06/23 13:37:39 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4Tubs.cc,v 1.74 2008/11/06 15:26:53 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    108108                "Invalid Z half-length");
    109109  }
    110   if ( pRMin < pRMax && pRMin >= 0 ) // Check radii
     110  if ( (pRMin < pRMax) && (pRMin >= 0) ) // Check radii
    111111  {
    112112    fRMin = pRMin ;
    … …  
    121121                "Invalid radii.");
    122122  }
    123   if ( pDPhi >= twopi ) // Check angles
     123
     124  fPhiFullTube = true;
     125  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
    124126  {
    125127    fDPhi=twopi;
     128    fSPhi=0;
    126129  }
    127130  else
    128131  {
     132    fPhiFullTube = false;
    129133    if ( pDPhi > 0 )
    130134    {
    … …  
    136140             << "        Negative delta-Phi ! - "
    137141             << pDPhi << G4endl;
    138       G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
    139                   "Invalid dphi.");
    140     }
    141   }
    142  
    143   // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
    144 
    145   fSPhi = pSPhi;
    146 
    147   if ( fSPhi < 0 )
    148   {
    149     fSPhi = twopi - std::fmod(std::fabs(fSPhi),twopi) ;
    150   }
    151   else
    152   {
    153     fSPhi = std::fmod(fSPhi,twopi) ;
    154   }
    155   if (fSPhi + fDPhi > twopi )
    156   {
    157     fSPhi -= twopi ;
    158   }
     142      G4Exception("G4Tubs::G4Tubs()", "InvalidSetup",
     143                  FatalException, "Invalid dphi.");
     144    }
     145
     146    // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
     147
     148    if ( pSPhi < 0 )
     149    {
     150      fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi);
     151    }
     152    else
     153    {
     154      fSPhi = std::fmod(pSPhi,twopi) ;
     155    }
     156    if ( fSPhi+fDPhi > twopi )
     157    {
     158      fSPhi -= twopi ;
     159    }
     160  }
     161  InitializeTrigonometry();
    159162}
    160163
    … …  
    290293        yoff2 = yMax    - yoffset ;
    291294
    292         if ( yoff1 >= 0 && yoff2 >= 0 ) // Y limits cross max/min x => no change
    293         {
     295        if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
     296        {                                   // => no change
    294297          pMin = xMin ;
    295298          pMax = xMax ;
    … …  
    317320        xoff2 = xMax - xoffset ;
    318321
    319         if ( xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change
    320         {
     322        if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
     323        {                                   // => no change
    321324          pMin = yMin ;
    322325          pMax = yMax ;
    … …  
    363366    noEntries = vertices->size() ;
    364367    noBetweenSections4 = noEntries - 4 ;
    365     /*
    366     G4cout << "vertices = " << noEntries << "\t"
    367            << "v-4 = " << noBetweenSections4 << G4endl;
    368     G4cout << G4endl;
    369     for(i = 0 ; i < noEntries ; i++ )
    370     {
    371       G4cout << i << "\t" << "v.x = " << ((*vertices)[i]).x() << "\t"
    372                           << "v.y = " << ((*vertices)[i]).y() << "\t"
    373                           << "v.z = " << ((*vertices)[i]).z() << "\t" << G4endl;
    374     }     
    375     G4cout << G4endl;
    376     G4cout << "ClipCrossSection" << G4endl;
    377     */
     368   
    378369    for (i = 0 ; i < noEntries ; i += 4 )
    379370    {
    380       // G4cout << "section = " << i << G4endl;
    381       ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
    382     }
    383     // G4cout << "ClipBetweenSections" << G4endl;
     371      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
     372    }
    384373    for (i = 0 ; i < noBetweenSections4 ; i += 4 )
    385374    {
    386       // G4cout << "between sections = " << i << G4endl;
    387       ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
    388     }
    389     if (pMin != kInfinity || pMax != -kInfinity )
     375      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
     376    }
     377    if ((pMin != kInfinity) || (pMax != -kInfinity) )
    390378    {
    391379      existsAfterClip = true ;
    … …  
    426414  G4double r2,pPhi,tolRMin,tolRMax;
    427415  EInside in = kOutside ;
    428  
    429   if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5)
     416  static const G4double halfCarTolerance=kCarTolerance*0.5;
     417  static const G4double halfRadTolerance=kRadTolerance*0.5;
     418  static const G4double halfAngTolerance=kAngTolerance*0.5;
     419
     420  if (std::fabs(p.z()) <= fDz - halfCarTolerance)
    430421  {
    431422    r2 = p.x()*p.x() + p.y()*p.y() ;
    432423
    433     if (fRMin) { tolRMin = fRMin + kRadTolerance*0.5 ; }
     424    if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
    434425    else       { tolRMin = 0 ; }
    435426
    436     tolRMax = fRMax - kRadTolerance*0.5 ;
     427    tolRMax = fRMax - halfRadTolerance ;
    437428     
    438     if (r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
    439     {
    440       if ( fDPhi == twopi )
     429    if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
     430    {
     431      if ( fPhiFullTube )
    441432      {
    442433        in = kInside ;
    … …  
    447438        // if not inside, try outer tolerant phi boundaries
    448439
    449         pPhi = std::atan2(p.y(),p.x()) ;
    450         if ((tolRMin==0)&&(p.x()==0)&&(p.y()==0))
     440        if ((tolRMin==0)&&(p.x()<=halfCarTolerance)&&(p.y()<=halfCarTolerance))
    451441        {
    452442          in=kSurface;
    … …  
    454444        else
    455445        {
    456           if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
     446          pPhi = std::atan2(p.y(),p.x()) ;
     447          if ( pPhi < -halfAngTolerance )  { pPhi += twopi; } // 0<=pPhi<2pi
    457448
    458449          if ( fSPhi >= 0 )
    459450          {
    460             if ( (std::abs(pPhi) < kAngTolerance*0.5)
    461               && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
     451            if ( (std::abs(pPhi) < halfAngTolerance)
     452              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
    462453            {
    463454              pPhi += twopi ; // 0 <= pPhi < 2pi
    464455            }
    465             if ( (pPhi >= fSPhi + kAngTolerance*0.5)
    466               && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) )
     456            if ( (pPhi >= fSPhi + halfAngTolerance)
     457              && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
    467458            {
    468459              in = kInside ;
    469460            }
    470             else if ( (pPhi >= fSPhi - kAngTolerance*0.5)
    471                    && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
     461            else if ( (pPhi >= fSPhi - halfAngTolerance)
     462                   && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
    472463            {
    473464              in = kSurface ;
    … …  
    476467          else  // fSPhi < 0
    477468          {
    478             if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
    479               && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
    480             else if ( (pPhi <= fSPhi + twopi + kAngTolerance*0.5)
    481                    && (pPhi >= fSPhi + fDPhi  - kAngTolerance*0.5) )
     469            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
     470              && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;} //kOutside
     471            else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
     472                   && (pPhi >= fSPhi + fDPhi  - halfAngTolerance) )
    482473            {
    483474              in = kSurface ;
    … …  
    493484    else  // Try generous boundaries
    494485    {
    495       tolRMin = fRMin - kRadTolerance*0.5 ;
    496       tolRMax = fRMax + kRadTolerance*0.5 ;
     486      tolRMin = fRMin - halfRadTolerance ;
     487      tolRMax = fRMax + halfRadTolerance ;
    497488
    498489      if ( tolRMin < 0 )  { tolRMin = 0; }
    … …  
    500491      if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
    501492      {
    502         if ( fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
    503         {
     493        if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
     494        {                        // Continuous in phi or on z-axis
    504495          in = kSurface ;
    505496        }
    … …  
    508499          pPhi = std::atan2(p.y(),p.x()) ;
    509500
    510           if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
     501          if ( pPhi < -halfAngTolerance)  { pPhi += twopi; } // 0<=pPhi<2pi
    511502          if ( fSPhi >= 0 )
    512503          {
    513             if ( (std::abs(pPhi) < kAngTolerance*0.5)
    514               && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
     504            if ( (std::abs(pPhi) < halfAngTolerance)
     505              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
    515506            {
    516507              pPhi += twopi ; // 0 <= pPhi < 2pi
    517508            }
    518             if ( (pPhi >= fSPhi - kAngTolerance*0.5)
    519               && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
     509            if ( (pPhi >= fSPhi - halfAngTolerance)
     510              && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
    520511            {
    521512              in = kSurface ;
    … …  
    524515          else  // fSPhi < 0
    525516          {
    526             if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
    527               && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
     517            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
     518              && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
    528519            else
    529520            {
    … …  
    535526    }
    536527  }
    537   else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5)
     528  else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
    538529  {                                          // Check within tolerant r limits
    539530    r2      = p.x()*p.x() + p.y()*p.y() ;
    540     tolRMin = fRMin - kRadTolerance*0.5 ;
    541     tolRMax = fRMax + kRadTolerance*0.5 ;
     531    tolRMin = fRMin - halfRadTolerance ;
     532    tolRMax = fRMax + halfRadTolerance ;
    542533
    543534    if ( tolRMin < 0 )  { tolRMin = 0; }
    … …  
    545536    if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
    546537    {
    547       if (fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
    548       {
     538      if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
     539      {                        // Continuous in phi or on z-axis
    549540        in = kSurface ;
    550541      }
    … …  
    553544        pPhi = std::atan2(p.y(),p.x()) ;
    554545
    555         if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; }  // 0<=pPhi<2pi
     546        if ( pPhi < -halfAngTolerance )  { pPhi += twopi; }  // 0<=pPhi<2pi
    556547        if ( fSPhi >= 0 )
    557548        {
    558           if ( (std::abs(pPhi) < kAngTolerance*0.5)
    559             && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
     549          if ( (std::abs(pPhi) < halfAngTolerance)
     550            && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
    560551          {
    561552            pPhi += twopi ; // 0 <= pPhi < 2pi
    562553          }
    563           if ( (pPhi >= fSPhi - kAngTolerance*0.5)
    564             && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
     554          if ( (pPhi >= fSPhi - halfAngTolerance)
     555            && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
    565556          {
    566557            in = kSurface;
    … …  
    569560        else  // fSPhi < 0
    570561        {
    571           if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
    572             && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
     562          if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
     563            && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;}
    573564          else
    574565          {
    … …  
    589580
    590581G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
    591 { G4int noSurfaces = 0;
     582{
     583  G4int noSurfaces = 0;
    592584  G4double rho, pPhi;
    593   G4double delta   = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
    594585  G4double distZ, distRMin, distRMax;
    595586  G4double distSPhi = kInfinity, distEPhi = kInfinity;
     587
     588  static const G4double halfCarTolerance = 0.5*kCarTolerance;
     589  static const G4double halfAngTolerance = 0.5*kAngTolerance;
     590
    596591  G4ThreeVector norm, sumnorm(0.,0.,0.);
    597592  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
    … …  
    604599  distZ    = std::fabs(std::fabs(p.z()) - fDz);
    605600
    606   if (fDPhi < twopi)   //  &&  rho ) // Protected against (0,0,z)
    607   {
    608     if ( rho )
     601  if (!fPhiFullTube)    // Protected against (0,0,z)
     602  {
     603    if ( rho > halfCarTolerance )
    609604    {
    610605      pPhi = std::atan2(p.y(),p.x());
    611606   
    612       if(pPhi  < fSPhi-delta)           { pPhi += twopi; }
    613       else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
    614 
    615       distSPhi = std::fabs( pPhi - fSPhi );       
     607      if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
     608      else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
     609
     610      distSPhi = std::fabs(pPhi - fSPhi);       
    616611      distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
    617612    }
    … …  
    624619    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
    625620  }
    626   if ( rho > delta )  { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
    627 
    628   if( distRMax <= delta )
     621  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
     622
     623  if( distRMax <= halfCarTolerance )
    629624  {
    630625    noSurfaces ++;
    631626    sumnorm += nR;
    632627  }
    633   if( fRMin && distRMin <= delta )
     628  if( fRMin && (distRMin <= halfCarTolerance) )
    634629  {
    635630    noSurfaces ++;
    … …  
    638633  if( fDPhi < twopi )   
    639634  {
    640     if (distSPhi <= dAngle)  // delta)
     635    if (distSPhi <= halfAngTolerance) 
    641636    {
    642637      noSurfaces ++;
    643638      sumnorm += nPs;
    644639    }
    645     if (distEPhi <= dAngle) // delta)
     640    if (distEPhi <= halfAngTolerance)
    646641    {
    647642      noSurfaces ++;
    … …  
    649644    }
    650645  }
    651   if (distZ <= delta) 
     646  if (distZ <= halfCarTolerance) 
    652647  {
    653648    noSurfaces ++;
    … …  
    668663  else if ( noSurfaces == 1 )  { norm = sumnorm; }
    669664  else                         { norm = sumnorm.unit(); }
     665
    670666  return norm;
    671667}
    … …  
    714710    }
    715711  }   
    716   if (fDPhi < twopi  &&  rho ) // Protected against (0,0,z)
     712  if (!fPhiFullTube  &&  rho ) // Protected against (0,0,z)
    717713  {
    718714    phi = std::atan2(p.y(),p.x()) ;
    … …  
    749745    case kNRMin : // Inner radius
    750746    {                     
    751       norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,0) ;
     747      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
    752748      break ;
    753749    }
    754750    case kNRMax : // Outer radius
    755751    {                 
    756       norm = G4ThreeVector(p.x()/rho,p.y()/rho,0) ;
     752      norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
    757753      break ;
    758754    }
    759755    case kNZ : //    + or - dz
    760756    {                             
    761       if ( p.z() > 0 ) norm = G4ThreeVector(0,0,1)  ;
    762       else             norm = G4ThreeVector(0,0,-1) ;
     757      if ( p.z() > 0 )  { norm = G4ThreeVector(0,0,1) ; }
     758      else              { norm = G4ThreeVector(0,0,-1); }
    763759      break ;
    764760    }
    765761    case kNSPhi:
    766762    {
    767       norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
     763      norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
    768764      break ;
    769765    }
    770766    case kNEPhi:
    771767    {
    772       norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
     768      norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
    773769      break;
    774770    }
    … …  
    804800//
    805801// NOTE:
    806 // - Precalculations for phi trigonometry are Done `just in time'
    807 // - `if valid' implies tolerant checking of intersection points
     802// - 'if valid' implies tolerant checking of intersection points
    808803
    809804G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
    810805                               const G4ThreeVector& v  ) const
    811806{
    812   G4double snxt = kInfinity ;  // snxt = default return value
    813 
    814   // Precalculated trig for phi intersections - used by r,z intersections to
    815   //                                            check validity
    816 
    817   G4bool seg ;        // true if segmented
    818 
    819   G4double hDPhi, hDPhiOT, hDPhiIT, cosHDPhiOT=0., cosHDPhiIT=0. ;
    820           // half dphi + outer tolerance
    821 
    822   G4double cPhi, sinCPhi=0., cosCPhi=0. ;  // central phi
    823 
    824   G4double tolORMin2, tolIRMax2 ;  // `generous' radii squared
    825 
     807  G4double snxt = kInfinity ;      // snxt = default return value
     808  G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
    826809  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
     810
     811  static const G4double halfCarTolerance = 0.5*kCarTolerance;
     812  static const G4double halfRadTolerance = 0.5*kRadTolerance;
    827813
    828814  // Intersection point variables
    829815  //
    830   G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi ;
    831 
    832   G4double t1, t2, t3, b, c, d ;   // Quadratic solver variables
    833 
    834   G4double Comp ;
    835   G4double cosSPhi, sinSPhi ;    // Trig for phi start intersect
    836 
    837   G4double ePhi, cosEPhi, sinEPhi ;  // for phi end intersect
    838 
    839   // Set phi divided flag and precalcs
    840 
    841   if ( fDPhi < twopi )
    842   {
    843     seg        = true ;
    844     hDPhi      = 0.5*fDPhi ;    // half delta phi
    845     cPhi       = fSPhi + hDPhi ;
    846     hDPhiOT    = hDPhi + 0.5*kAngTolerance ;  // outers tol' half delta phi
    847     hDPhiIT    = hDPhi - 0.5*kAngTolerance ;
    848     sinCPhi    = std::sin(cPhi) ;
    849     cosCPhi    = std::cos(cPhi) ;
    850     cosHDPhiOT = std::cos(hDPhiOT) ;
    851     cosHDPhiIT = std::cos(hDPhiIT) ;
    852   }
    853   else
    854   {
    855     seg = false  ;
    856   }
    857 
     816  G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
     817  G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables
     818 
    858819  // Calculate tolerant rmin and rmax
    859820
    860821  if (fRMin > kRadTolerance)
    861822  {
    862     tolORMin2 = (fRMin - 0.5*kRadTolerance)*(fRMin - 0.5*kRadTolerance) ;
    863     tolIRMin2 = (fRMin + 0.5*kRadTolerance)*(fRMin + 0.5*kRadTolerance) ;
     823    tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
     824    tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
    864825  }
    865826  else
    … …  
    868829    tolIRMin2 = 0.0 ;
    869830  }
    870   tolORMax2 = (fRMax + 0.5*kRadTolerance)*(fRMax + 0.5*kRadTolerance) ;
    871   tolIRMax2 = (fRMax - 0.5*kRadTolerance)*(fRMax - 0.5*kRadTolerance) ;
     831  tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
     832  tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
    872833
    873834  // Intersection with Z surfaces
    874835
    875   tolIDz = fDz - kCarTolerance*0.5 ;
    876   tolODz = fDz + kCarTolerance*0.5 ;
     836  tolIDz = fDz - halfCarTolerance ;
     837  tolODz = fDz + halfCarTolerance ;
    877838
    878839  if (std::fabs(p.z()) >= tolIDz)
    … …  
    880841    if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
    881842    {
    882       s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;     // Z intersect distance
     843      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;   // Z intersect distance
    883844
    884845      if(s < 0.0)  { s = 0.0; }
    … …  
    890851      // Check validity of intersection
    891852
    892       if (tolIRMin2 <= rho2 && rho2 <= tolIRMax2)
    893       {
    894         if (seg && rho2)
     853      if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
     854      {
     855        if (!fPhiFullTube && rho2)
    895856        {
    896857          // Psi = angle made with central (average) phi of shape
    … …  
    899860          iden   = std::sqrt(rho2) ;
    900861          cosPsi = inum/iden ;
    901           if (cosPsi >= cosHDPhiIT) return s ;
     862          if (cosPsi >= cosHDPhiIT)  { return s ; }
    902863        }
    903864        else
    … …  
    909870    else
    910871    {
    911       if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
     872      if ( snxt<halfCarTolerance )  { snxt=0; }
    912873      return snxt ;  // On/outside extent, and heading away
    913874                     // -> cannot intersect
    … …  
    934895    b = t2/t1 ;
    935896    c = t3 - fRMax*fRMax ;
    936     if (t3 >= tolORMax2 && t2<0)   // This also handles the tangent case
     897    if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
    937898    {
    938899      // Try outer cylinder intersection
    … …  
    954915            // Z ok. Check phi intersection if reqd
    955916            //
    956             if (!seg)
     917            if (fPhiFullTube)
    957918            {
    958919              return s ;
    … …  
    963924              yi     = p.y() + s*v.y() ;
    964925              cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
    965               if (cosPsi >= cosHDPhiIT) return s ;
     926              if (cosPsi >= cosHDPhiIT)  { return s ; }
    966927            }
    967928          }  //  end if std::fabs(zi)
    … …  
    974935      // check not inside, and heading through tubs (-> 0 to in)
    975936
    976       if (t3 > tolIRMin2 && t2 < 0 && std::fabs(p.z()) <= tolIDz)
     937      if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
    977938      {
    978939        // Inside both radii, delta r -ve, inside z extent
    979940
    980         if (seg)
     941        if (!fPhiFullTube)
    981942        {
    982943          inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
    … …  
    1004965                snxt = c/(-b+std::sqrt(d)); // using safe solution
    1005966                                            // for quadratic equation
    1006                 if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
     967                if ( snxt < halfCarTolerance ) { snxt=0; }
    1007968                return snxt ;
    1008969              }     
    … …  
    1035996              snxt= c/(-b+std::sqrt(d)); // using safe solution
    1036997                                         // for quadratic equation
    1037               if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
     998              if ( snxt < halfCarTolerance ) { snxt=0; }
    1038999              return snxt ;
    10391000            }     
    … …  
    10431004            }
    10441005          }
    1045         } // end if   (seg)
     1006        } // end if   (!fPhiFullTube)
    10461007      }   // end if   (t3>tolIRMin2)
    10471008    }     // end if   (Inside Outer Radius)
    … …  
    10561017
    10571018        s = -b + std::sqrt(d) ;
    1058         if (s >= -0.5*kCarTolerance)  // check forwards
     1019        if (s >= -halfCarTolerance)  // check forwards
    10591020        {
    10601021          // Check z intersection
    … …  
    10661027            // Z ok. Check phi
    10671028            //
    1068             if ( !seg )
     1029            if ( fPhiFullTube )
    10691030            {
    10701031              return s ;
    … …  
    10981059  //         -> use some form of loop Construct ?
    10991060  //
    1100   if ( seg )
     1061  if ( !fPhiFullTube )
    11011062  {
    11021063    // First phi surface (Starting phi)
    1103 
    1104     sinSPhi = std::sin(fSPhi) ;
    1105     cosSPhi = std::cos(fSPhi) ;
     1064    //
    11061065    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
    11071066                   
    … …  
    11101069      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
    11111070
    1112       if ( Dist < kCarTolerance*0.5 )
     1071      if ( Dist < halfCarTolerance )
    11131072      {
    11141073        s = Dist/Comp ;
    … …  
    11351094              // - check intersecting with correct half-plane
    11361095              //
    1137               if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;
    1138             }   
     1096              if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = s; }
     1097            }
    11391098          }
    11401099        }
    … …  
    11421101    }
    11431102     
    1144     // Second phi surface (`E'nding phi)
    1145 
    1146     ePhi    = fSPhi + fDPhi ;
    1147     sinEPhi = std::sin(ePhi) ;
    1148     cosEPhi = std::cos(ePhi) ;
     1103    // Second phi surface (Ending phi)
     1104
    11491105    Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
    11501106       
    … …  
    11531109      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
    11541110
    1155       if ( Dist < kCarTolerance*0.5 )
     1111      if ( Dist < halfCarTolerance )
    11561112      {
    11571113        s = Dist/Comp ;
    … …  
    11771133              // - check intersecting with correct half-plane
    11781134              //
    1179               if ( (yi*cosCPhi-xi*sinCPhi) >= 0 )  { snxt = s; }
    1180             }   
     1135              if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = s; }
     1136            }                         //?? >=-halfCarTolerance
    11811137          }
    11821138        }
    11831139      }
    11841140    }         //  Comp < 0
    1185   }           //  seg != 0
    1186   if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
     1141  }           //  !fPhiFullTube
     1142  if ( snxt<halfCarTolerance )  { snxt=0; }
    11871143  return snxt ;
    11881144}
    … …  
    12171173{
    12181174  G4double safe=0.0, rho, safe1, safe2, safe3 ;
    1219   G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
     1175  G4double safePhi, cosPsi ;
    12201176
    12211177  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
    … …  
    12281184  if ( safe3 > safe )  { safe = safe3; }
    12291185
    1230   if (fDPhi < twopi && rho)
    1231   {
    1232     phiC    = fSPhi + fDPhi*0.5 ;
    1233     cosPhiC = std::cos(phiC) ;
    1234     sinPhiC = std::sin(phiC) ;
    1235 
     1186  if ( (!fPhiFullTube) && (rho) )
     1187  {
    12361188    // Psi=angle from central phi to point
    12371189    //
    1238     cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
    1239 
     1190    cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
     1191   
    12401192    if ( cosPsi < std::cos(fDPhi*0.5) )
    12411193    {
    12421194      // Point lies outside phi range
    12431195
    1244       if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
    1245       {
    1246         safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
     1196      if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
     1197      {
     1198        safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
    12471199      }
    12481200      else
    12491201      {
    1250         ePhi    = fSPhi + fDPhi ;
    1251         safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
     1202        safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
    12521203      }
    12531204      if ( safePhi > safe )  { safe = safePhi; }
    … …  
    12691220                                      G4ThreeVector *n    ) const
    12701221{ 
    1271   ESide side = kNull , sider = kNull, sidephi = kNull ;
    1272   G4double snxt, sr = kInfinity, sphi = kInfinity, pdist ;
     1222  ESide side=kNull , sider=kNull, sidephi=kNull ;
     1223  G4double snxt, sr=kInfinity, sphi=kInfinity, pdist ;
    12731224  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
    12741225
     1226  static const G4double halfCarTolerance = kCarTolerance*0.5;
     1227  static const G4double halfAngTolerance = kAngTolerance*0.5;
     1228 
    12751229  // Vars for phi intersection:
    12761230
    1277   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
    1278   G4double cPhi, sinCPhi, cosCPhi ;
    12791231  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
    12801232 
    … …  
    12841236  {
    12851237    pdist = fDz - p.z() ;
    1286     if ( pdist > kCarTolerance*0.5 )
     1238    if ( pdist > halfCarTolerance )
    12871239    {
    12881240      snxt = pdist/v.z() ;
    … …  
    13031255    pdist = fDz + p.z() ;
    13041256
    1305     if ( pdist > kCarTolerance*0.5 )
     1257    if ( pdist > halfCarTolerance )
    13061258    {
    13071259      snxt = -pdist/v.z() ;
    … …  
    13261278  // Radial Intersections
    13271279  //
    1328   // Find intersction with cylinders at rmax/rmin
     1280  // Find intersection with cylinders at rmax/rmin
    13291281  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
    13301282  //
    … …  
    13591311        b     = t2/t1 ;
    13601312        c     = deltaR/t1 ;
    1361         d2= b*b-c;
    1362         if(d2>=0.){sr    = -b + std::sqrt(d2);}
    1363         else{sr=0.;};
     1313        d2    = b*b-c;
     1314        if( d2 >= 0 ) { sr = -b + std::sqrt(d2); }
     1315        else          { sr = 0.; }
    13641316        sider = kRMax ;
    13651317      }
    … …  
    13711323        if ( calcNorm )
    13721324        {
    1373           // if ( p.x() || p.y() )
    1374           // {
    1375           //  *n=G4ThreeVector(p.x(),p.y(),0);
    1376           // }
    1377           // else
    1378           // {
    1379           //  *n=v;
    1380           // }
    13811325          *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
    13821326          *validNorm = true ;
    … …  
    14171361          c     = deltaR/t1 ;
    14181362          d2    = b*b-c;
    1419           if(d2>=0.)
     1363          if( d2 >=0. )
    14201364          {
    14211365            sr     = -b + std::sqrt(d2) ;
    … …  
    14231367          }
    14241368          else // Case: On the border+t2<kRadTolerance
    1425                //       (v is perpendiculair to the surface)
     1369               //       (v is perpendicular to the surface)
    14261370          {
    14271371            if (calcNorm)
    … …  
    14411385        c      = deltaR/t1;
    14421386        d2     = b*b-c;
    1443         if(d2>=0.)
     1387        if( d2 >= 0 )
    14441388        {
    14451389          sr     = -b + std::sqrt(d2) ;
    … …  
    14471391        }
    14481392        else // Case: On the border+t2<kRadTolerance
    1449              //       (v is perpendiculair to the surface)
     1393             //       (v is perpendicular to the surface)
    14501394        {
    14511395          if (calcNorm)
    … …  
    14611405    // Phi Intersection
    14621406
    1463     if ( fDPhi < twopi )
    1464     {
    1465       sinSPhi = std::sin(fSPhi) ;
    1466       cosSPhi = std::cos(fSPhi) ;
    1467       ePhi    = fSPhi + fDPhi ;
    1468       sinEPhi = std::sin(ePhi) ;
    1469       cosEPhi = std::cos(ePhi) ;
    1470       cPhi    = fSPhi + fDPhi*0.5 ;
    1471       sinCPhi = std::sin(cPhi) ;
    1472       cosCPhi = std::cos(cPhi) ;
    1473 
     1407    if ( !fPhiFullTube )
     1408    {
    14741409      // add angle calculation with correction
    14751410      // of the difference in domain of atan2 and Sphi
    14761411      //
    14771412      vphi = std::atan2(v.y(),v.x()) ;
    1478 
    1479       if ( vphi < fSPhi - kAngTolerance*0.5  )             { vphi += twopi; }
    1480       else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) { vphi -= twopi; }
     1413     
     1414      if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
     1415      else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
    14811416
    14821417
    … …  
    14921427        compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
    14931428        compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
     1429       
    14941430        sidephi = kNull;
    14951431       
    1496         if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
    1497                               && (pDistE <= 0.5*kCarTolerance) ) )
    1498          || ( (fDPhi >  pi) && !((pDistS >  0.5*kCarTolerance)
    1499                               && (pDistE >  0.5*kCarTolerance) ) )  )
     1432        if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
     1433                              && (pDistE <= halfCarTolerance) ) )
     1434         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
     1435                              && (pDistE >  halfCarTolerance) ) )  )
    15001436        {
    15011437          // Inside both phi *full* planes
    … …  
    15051441            sphi = pDistS/compS ;
    15061442           
    1507             if (sphi >= -0.5*kCarTolerance)
     1443            if (sphi >= -halfCarTolerance)
    15081444            {
    15091445              xi = p.x() + sphi*v.x() ;
    … …  
    15131449              // (if not -> no intersect)
    15141450              //
    1515               if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
    1516                 { sidephi = kSPhi;
    1517                 if (((fSPhi-0.5*kAngTolerance)<=vphi)
    1518                    &&((ePhi+0.5*kAngTolerance)>=vphi))
     1451              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
     1452              {
     1453                sidephi = kSPhi;
     1454                if (((fSPhi-halfAngTolerance)<=vphi)
     1455                   &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
    15191456                {
    15201457                  sphi = kInfinity;
    15211458                }
    15221459              }
    1523               else if ((yi*cosCPhi-xi*sinCPhi)>=0)
     1460              else if ( yi*cosCPhi-xi*sinCPhi >=0 )
    15241461              {
    15251462                sphi = kInfinity ;
    … …  
    15281465              {
    15291466                sidephi = kSPhi ;
    1530                 if ( pDistS > -kCarTolerance*0.5 )
     1467                if ( pDistS > -halfCarTolerance )
    15311468                {
    15321469                  sphi = 0.0 ; // Leave by sphi immediately
    … …  
    15501487            // Only check further if < starting phi intersection
    15511488            //
    1552             if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
     1489            if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
    15531490            {
    15541491              xi = p.x() + sphi2*v.x() ;
    … …  
    15591496                // Leaving via ending phi
    15601497                //
    1561                 if(!(((fSPhi-0.5*kAngTolerance)<=vphi)
    1562                     &&((ePhi+0.5*kAngTolerance)>=vphi)))
     1498                if( !((fSPhi-halfAngTolerance <= vphi)
     1499                     &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
    15631500                {
    15641501                  sidephi = kEPhi ;
    1565                   if ( pDistE <= -kCarTolerance*0.5 )  { sphi = sphi2 ; }
    1566                   else                                 { sphi = 0.0 ; }
     1502                  if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
     1503                  else                                { sphi = 0.0 ;  }
    15671504                }
    15681505              }
    … …  
    15741511                //
    15751512                sidephi = kEPhi ;
    1576                 if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
    1577                 else                                { sphi = 0.0 ; }
     1513                if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
     1514                else                               { sphi = 0.0 ;  }
    15781515              }
    15791516            }
    … …  
    15891526        // On z axis + travel not || to z axis -> if phi of vector direction
    15901527        // within phi of shape, Step limited by rmax, else Step =0
    1591 
    1592         // vphi = std::atan2(v.y(),v.x()) ;//defined previosly
    1593         // G4cout<<"In axis vphi="<<vphi
    1594         //       <<" Sphi="<<fSPhi<<" Ephi="<<ePhi<<G4endl;
    1595         // old  if ( (fSPhi < vphi) && (vphi < fSPhi + fDPhi) )
    1596         // new : correction for if statement, must be '<=' 
    15971528               
    1598         if ( ((fSPhi-0.5*kAngTolerance) <= vphi)
    1599            && (vphi <=( ePhi+0.5*kAngTolerance) ))
     1529        if ( (fSPhi - halfAngTolerance <= vphi)
     1530           && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
    16001531        {
    16011532          sphi = kInfinity ;
    … …  
    16401571        if ( fDPhi <= pi )
    16411572        {
    1642           *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
     1573          *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
    16431574          *validNorm = true ;
    16441575        }
    … …  
    16521583        if (fDPhi <= pi)
    16531584        {
    1654           *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
     1585          *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
    16551586          *validNorm = true ;
    16561587        }
    … …  
    16621593
    16631594      case kPZ:
    1664         *n=G4ThreeVector(0,0,1) ;
    1665         *validNorm=true ;
     1595        *n         = G4ThreeVector(0,0,1) ;
     1596        *validNorm = true ;
    16661597        break ;
    16671598
    … …  
    16901621    }
    16911622  }
    1692   if ( snxt<kCarTolerance*0.5 )  { snxt=0 ; }
     1623  if ( snxt<halfCarTolerance )  { snxt=0 ; }
    16931624
    16941625  return snxt ;
    … …  
    17011632G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
    17021633{
    1703   G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
    1704   G4double safePhi, phiC, cosPhiC, sinPhiC, ePhi ;
     1634  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
    17051635  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
    17061636
    … …  
    17381668  // Check if phi divided, Calc distances closest phi plane
    17391669  //
    1740   if ( fDPhi < twopi )
    1741   {
    1742     // Above/below central phi of Tubs?
    1743 
    1744     phiC    = fSPhi + fDPhi*0.5 ;
    1745     cosPhiC = std::cos(phiC) ;
    1746     sinPhiC = std::sin(phiC) ;
    1747 
    1748     if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
    1749     {
    1750       safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
     1670  if ( !fPhiFullTube )
     1671  {
     1672    if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
     1673    {
     1674      safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
    17511675    }
    17521676    else
    17531677    {
    1754       ePhi    = fSPhi + fDPhi ;
    1755       safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
     1678      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
    17561679    }
    17571680    if (safePhi < safe)  { safe = safePhi ; }
    … …  
    18061729  // on the x axis. Will give better extent calculations when not rotated.
    18071730
    1808   if (fDPhi == pi*2.0 && fSPhi == 0 )  { sAngle = -meshAngle*0.5 ; }
    1809   else                                 { sAngle =  fSPhi ; }
     1731  if (fPhiFullTube && (fSPhi == 0) )  { sAngle = -meshAngle*0.5 ; }
     1732  else                                { sAngle =  fSPhi ; }
    18101733   
    18111734  vertices = new G4ThreeVectorList();
    … …  
    19751898  if (fRMin != 0)
    19761899  {
    1977     if (fDPhi >= twopi)
     1900    if (fPhiFullTube)
    19781901    {
    19791902      pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
    … …  
    19861909  else
    19871910  {
    1988     if (fDPhi >= twopi)
     1911    if (fPhiFullTube)
    19891912    {
    19901913      pNURBS = new G4NURBScylinder (fRMax,fDz) ;
Note: See TracChangeset for help on using the changeset viewer.