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/G4Cons.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4Cons.cc,v 1.56 2008/02/20 08:56:16 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4Cons.cc,v 1.60 2008/11/06 15:26:53 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    8787
    8888  if ( pDz > 0 )
    89     fDz = pDz ;
     89  {
     90    fDz = pDz;
     91  }
    9092  else
    9193  {
    … …  
    99101  // Check radii
    100102
    101   if ( pRmin1 < pRmax1 && pRmin2 < pRmax2 && pRmin1 >= 0 && pRmin2 >= 0 )
     103  if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) )
    102104  {
    103105
    … …  
    106108    fRmin2 = pRmin2 ;
    107109    fRmax2 = pRmax2 ;
    108     if( (pRmin1 == 0.0 && pRmin2 > 0.0) ) fRmin1 = 1e3*kRadTolerance ;
    109     if( (pRmin2 == 0.0 && pRmin1 > 0.0) ) fRmin2 = 1e3*kRadTolerance ;
     110    if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
     111    if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
    110112  }
    111113  else
    … …  
    119121  }
    120122
    121   // Check angles
    122 
    123   if ( pDPhi >= twopi )
     123  fPhiFullCone = true;
     124  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
    124125  {
    125126    fDPhi=twopi;
    … …  
    128129  else
    129130  {
    130     if ( pDPhi > 0 ) fDPhi = pDPhi ;
     131    fPhiFullCone = false;
     132    if ( pDPhi > 0 )
     133    {
     134      fDPhi = pDPhi;
     135    }
    131136    else
    132137    {
    … …  
    135140             << pDPhi << G4endl;
    136141      G4Exception("G4Cons::G4Cons()", "InvalidSetup",
    137                   FatalException, "Invalid pDPhi.") ;
    138     }
    139 
    140     // Ensure pSPhi in 0-2PI or -2PI-0 range if shape crosses 0
    141 
    142     if ( pSPhi < 0 ) fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi) ;
    143     else             fSPhi = std::fmod(pSPhi,twopi) ;
    144      
    145     if (fSPhi + fDPhi > twopi) fSPhi -= twopi ;
    146   }
     142                  FatalException, "Invalid dphi.");
     143    }
     144
     145    // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
     146
     147    if ( pSPhi < 0 )
     148    {
     149      fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi);
     150    }
     151    else
     152    {
     153      fSPhi = std::fmod(pSPhi,twopi) ;
     154    }
     155    if ( fSPhi+fDPhi > twopi )
     156    {
     157      fSPhi -= twopi ;
     158    }
     159  }
     160  InitializeTrigonometry();
    147161}
    148162
    … …  
    173187  G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
    174188  EInside in;
    175 
    176   if (std::fabs(p.z()) > fDz + kCarTolerance*0.5 ) return in = kOutside;
    177   else if(std::fabs(p.z()) >= fDz - kCarTolerance*0.5 )   in = kSurface;
    178   else                                                    in = kInside;
     189  static const G4double halfCarTolerance=kCarTolerance*0.5;
     190  static const G4double halfRadTolerance=kRadTolerance*0.5;
     191  static const G4double halfAngTolerance=kAngTolerance*0.5;
     192
     193  if (std::fabs(p.z()) > fDz + halfCarTolerance )  { return in = kOutside; }
     194  else if(std::fabs(p.z()) >= fDz - halfCarTolerance )    { in = kSurface; }
     195  else                                                    { in = kInside;  }
    179196
    180197  r2 = p.x()*p.x() + p.y()*p.y() ;
    … …  
    184201  // rh2 = rh*rh;
    185202
    186   tolRMin = rl - kRadTolerance*0.5 ;
    187   if ( tolRMin < 0 ) tolRMin = 0 ;
    188   tolRMax = rh + kRadTolerance*0.5 ;
    189 
    190   if ( r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax ) return in = kOutside;
    191 
    192   if (rl) tolRMin = rl + kRadTolerance*0.5 ;
    193   else    tolRMin = 0.0 ;
    194   tolRMax         = rh - kRadTolerance*0.5 ;
     203  tolRMin = rl - halfRadTolerance;
     204  if ( tolRMin < 0 )  { tolRMin = 0; }
     205  tolRMax = rh + halfRadTolerance;
     206
     207  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
     208
     209  if (rl) { tolRMin = rl + halfRadTolerance; }
     210  else    { tolRMin = 0.0; }
     211  tolRMax = rh - halfRadTolerance;
    195212     
    196213  if (in == kInside) // else it's kSurface already
    197214  {
    198      if (r2 < tolRMin*tolRMin || r2 >= tolRMax*tolRMax) in = kSurface;
    199     //  if (r2 <= tolRMin*tolRMin || r2-rh2 >= -rh*kRadTolerance) in = kSurface;
    200   }
    201   if ( ( fDPhi < twopi - kAngTolerance ) &&
    202        ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
     215     if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
     216  }
     217  if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
    203218  {
    204219    pPhi = std::atan2(p.y(),p.x()) ;
    205220
    206     if ( pPhi < fSPhi - kAngTolerance*0.5  )             pPhi += twopi ;
    207     else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 ) pPhi -= twopi;
     221    if ( pPhi < fSPhi - halfAngTolerance  )             { pPhi += twopi; }
     222    else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
    208223   
    209     if ( (pPhi < fSPhi - kAngTolerance*0.5) ||         
    210          (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) )   return in = kOutside;
     224    if ( (pPhi < fSPhi - halfAngTolerance) ||         
     225         (pPhi > fSPhi + fDPhi + halfAngTolerance) )  { return in = kOutside; }
    211226     
    212227    else if (in == kInside)  // else it's kSurface anyway already
    213228    {
    214         if ( (pPhi < fSPhi + kAngTolerance*0.5) ||
    215              (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) )  in = kSurface ;
    216     }
    217   }
    218   else if( fDPhi < twopi - kAngTolerance )  in = kSurface ;
     229       if ( (pPhi < fSPhi + halfAngTolerance) ||
     230            (pPhi > fSPhi + fDPhi - halfAngTolerance) )  { in = kSurface; }
     231    }
     232  }
     233  else if ( !fPhiFullCone )  { in = kSurface; }
    219234
    220235  return in ;
    … …  
    244259                    G4double&          pMax  ) const
    245260{
    246   if ( !pTransform.IsRotated() &&
    247         fDPhi == twopi && fRmin1 == 0 && fRmin2 == 0 )
     261  if ( !pTransform.IsRotated() && (fDPhi == twopi)
     262    && (fRmin1 == 0) && (fRmin2 == 0) )
    248263  {
    249264    // Special case handling for unrotated solid cones
    … …  
    265280    if (pVoxelLimit.IsZLimited())
    266281    {
    267       if( zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance ||
    268           zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance   )
     282      if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
     283          (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)  )
    269284      {
    270285        return false ;
    … …  
    290305    if (pVoxelLimit.IsXLimited())
    291306    {
    292       if ( xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance ||
    293            xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance    )
     307      if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
     308           (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance)  )
    294309      {
    295310        return false ;
    … …  
    315330    if (pVoxelLimit.IsYLimited())
    316331    {
    317       if ( yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance ||
    318            yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance     )
     332      if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
     333           (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance)  )
    319334      {
    320335        return false ;
    … …  
    338353        yoff2 = yMax - yoffset ;
    339354
    340         if (yoff1 >= 0 && yoff2 >= 0) // Y limits cross max/min x => no change
    341         {
     355        if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
     356        {                                 // => no change
    342357          pMin = xMin ;
    343358          pMax = xMax ;
    … …  
    362377        xoff2 = xMax - xoffset ;
    363378
    364         if (xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change
    365         {
     379        if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
     380        {                                  // => no change
    366381          pMin = yMin ;
    367382          pMax = yMax ;
    … …  
    409424    for ( i = 0 ; i < noEntries ; i += 4 )
    410425    {
    411       ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
     426      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    412427    }
    413428    for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
    414429    {
    415       ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
     430      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    416431    }   
    417     if ( pMin != kInfinity || pMax != -kInfinity )
     432    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
    418433    {
    419434      existsAfterClip = true ;
    … …  
    462477  G4double tanRMin, secRMin, pRMin, widRMin;
    463478  G4double tanRMax, secRMax, pRMax, widRMax;
    464   G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
     479
     480  static const G4double delta  = 0.5*kCarTolerance;
     481  static const G4double dAngle = 0.5*kAngTolerance;
    465482 
    466   G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.0);
     483  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
    467484  G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
    468485
    … …  
    482499  distRMax = std::fabs(pRMax - widRMax)/secRMax;
    483500
    484   if (fDPhi < twopi)   //  &&  rho ) // Protected against (0,0,z)
     501  if (!fPhiFullCone)  // Protected against (0,0,z)
    485502  {
    486503    if ( rho )
    … …  
    488505      pPhi = std::atan2(p.y(),p.x());
    489506
    490       if(pPhi  < fSPhi-delta)           pPhi     += twopi;
    491       else if(pPhi > fSPhi+fDPhi+delta) pPhi     -= twopi;
     507      if (pPhi  < fSPhi-delta)            { pPhi += twopi; }
     508      else if (pPhi > fSPhi+fDPhi+delta)  { pPhi -= twopi; }
    492509
    493510      distSPhi = std::fabs( pPhi - fSPhi );
    494       distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
     511      distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
    495512    }
    496513    else if( !(fRmin1) || !(fRmin2) )
    … …  
    499516      distEPhi = 0.;
    500517    }
    501     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
    502     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
     518    nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
     519    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
    503520  }
    504521  if ( rho > delta )   
    505522  {
    506     nR = G4ThreeVector(p.x()/rho/secRMax,p.y()/rho/secRMax,-tanRMax/secRMax);
    507     if (fRmin1 || fRmin2) nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
     523    nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
     524    if (fRmin1 || fRmin2)
     525    {
     526      nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
     527    }
    508528  }
    509529
    … …  
    513533    sumnorm += nR;
    514534  }
    515   if( (fRmin1 || fRmin2) && distRMin <= delta )
     535  if( (fRmin1 || fRmin2) && (distRMin <= delta) )
    516536  {
    517537    noSurfaces ++;
    518538    sumnorm += nr;
    519539  }
    520   if( fDPhi < twopi )   
     540  if( !fPhiFullCone )   
    521541  {
    522542    if (distSPhi <= dAngle)
    … …  
    534554  {
    535555    noSurfaces ++;
    536     if ( p.z() >= 0.)  sumnorm += nZ;
    537     else               sumnorm -= nZ;
     556    if ( p.z() >= 0.)  { sumnorm += nZ; }
     557    else               { sumnorm -= nZ; }
    538558  }
    539559  if ( noSurfaces == 0 )
    … …  
    545565     norm = ApproxSurfaceNormal(p);
    546566  }
    547   else if ( noSurfaces == 1 ) norm = sumnorm;
    548   else                        norm = sumnorm.unit();
     567  else if ( noSurfaces == 1 )  { norm = sumnorm; }
     568  else                         { norm = sumnorm.unit(); }
     569
    549570  return norm ;
    550571}
    551572
    552 //////////////////////////////////////////////////////////////////////////////////
     573////////////////////////////////////////////////////////////////////////////
    553574//
    554575// Algorithm for SurfaceNormal() following the original specification
    … …  
    605626    }
    606627  }
    607   if ( fDPhi < twopi && rho )  // Protected against (0,0,z)
     628  if ( !fPhiFullCone && rho )  // Protected against (0,0,z)
    608629  {
    609630    phi = std::atan2(p.y(),p.x()) ;
    610631
    611     if (phi < 0) phi += twopi ;
    612 
    613     if (fSPhi < 0) distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
    614     else           distSPhi = std::fabs(phi - fSPhi)*rho ;
     632    if (phi < 0)  { phi += twopi; }
     633
     634    if (fSPhi < 0)  { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
     635    else            { distSPhi = std::fabs(phi - fSPhi)*rho; }
    615636
    616637    distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
    … …  
    620641    if (distSPhi < distEPhi)
    621642    {
    622       if (distSPhi < distMin) side = kNSPhi ;
     643      if (distSPhi < distMin)  { side = kNSPhi; }
    623644    }
    624645    else
    625646    {
    626       if (distEPhi < distMin) side = kNEPhi ;
     647      if (distEPhi < distMin)  { side = kNEPhi; }
    627648    }
    628649  }   
    … …  
    631652    case kNRMin:      // Inner radius
    632653      rho *= secRMin ;
    633       norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,tanRMin/secRMin) ;
     654      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
    634655      break ;
    635656    case kNRMax:      // Outer radius
    636657      rho *= secRMax ;
    637       norm = G4ThreeVector(p.x()/rho,p.y()/rho,-tanRMax/secRMax) ;
     658      norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
    638659      break ;
    639660    case kNZ:      // +/- dz
    640       if (p.z() > 0) norm = G4ThreeVector(0,0,1) ;
    641       else           norm = G4ThreeVector(0,0,-1) ;
     661      if (p.z() > 0)  { norm = G4ThreeVector(0,0,1);  }
     662      else            { norm = G4ThreeVector(0,0,-1); }
    642663      break ;
    643664    case kNSPhi:
    644       norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
     665      norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
    645666      break ;
    646667    case kNEPhi:
    647       norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
     668      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
    648669      break ;
    649670    default:
    … …  
    677698//
    678699// NOTE:
    679 // - Precalculations for phi trigonometry are Done `just in time'
    680700// - `if valid' implies tolerant checking of intersection points
    681701// - z, phi intersection from Tubs
    … …  
    686706  G4double snxt = kInfinity ;      // snxt = default return value
    687707
    688   G4bool seg ;        // true if segmented in phi
    689   G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0. ;
    690           // half dphi + outer tolerance
    691   G4double cPhi,sinCPhi=0.,cosCPhi=0. ;  // central phi
     708  static const G4double halfCarTolerance=kCarTolerance*0.5;
     709  static const G4double halfRadTolerance=kRadTolerance*0.5;
    692710
    693711  G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;  // Data for cones
    … …  
    704722  G4double nt1,nt2,nt3 ;
    705723  G4double Comp ;
    706   G4double cosSPhi,sinSPhi ;    // Trig for phi start intersect
    707   G4double ePhi,cosEPhi,sinEPhi ;  // for phi end intersect
    708 
    709   //
    710   // Set phi divided flag and precalcs
    711   //
    712   if (fDPhi < twopi)
    713   {
    714     seg        = true ;
    715     hDPhi      = 0.5*fDPhi ;    // half delta phi
    716     cPhi       = fSPhi + hDPhi ; ;
    717     hDPhiOT    = hDPhi + 0.5*kAngTolerance ;  // outers tol' half delta phi
    718     hDPhiIT    = hDPhi - 0.5*kAngTolerance ;
    719     sinCPhi    = std::sin(cPhi) ;
    720     cosCPhi    = std::cos(cPhi) ;
    721     cosHDPhiOT = std::cos(hDPhiOT) ;
    722     cosHDPhiIT = std::cos(hDPhiIT) ;
    723   }
    724   else     seg = false ;
    725724
    726725  // Cone Precalcs
    … …  
    730729  rMinAv  = (fRmin1 + fRmin2)*0.5 ;
    731730
    732   if (rMinAv > kRadTolerance*0.5)
    733   {
    734     rMinOAv = rMinAv - kRadTolerance*0.5 ;
    735     rMinIAv = rMinAv + kRadTolerance*0.5 ;
     731  if (rMinAv > halfRadTolerance)
     732  {
     733    rMinOAv = rMinAv - halfRadTolerance ;
     734    rMinIAv = rMinAv + halfRadTolerance ;
    736735  }
    737736  else
    … …  
    743742  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
    744743  rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
    745   rMaxOAv = rMaxAv + kRadTolerance*0.5 ;
     744  rMaxOAv = rMaxAv + halfRadTolerance ;
    746745   
    747746  // Intersection with z-surfaces
    748747
    749   tolIDz = fDz - kCarTolerance*0.5 ;
    750   tolODz = fDz + kCarTolerance*0.5 ;
     748  tolIDz = fDz - halfCarTolerance ;
     749  tolODz = fDz + halfCarTolerance ;
    751750
    752751  if (std::fabs(p.z()) >= tolIDz)
    … …  
    756755      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;  // Z intersect distance
    757756
    758       if( s < 0.0 ) s = 0.0 ;                  // negative dist -> zero
     757      if( s < 0.0 )  { s = 0.0; }                      // negative dist -> zero
    759758
    760759      xi   = p.x() + s*v.x() ;  // Intersection coords
    761760      yi   = p.y() + s*v.y() ;
    762       rhoi2 = xi*xi + yi*yi   ;
     761      rhoi2 = xi*xi + yi*yi  ;
    763762
    764763      // Check validity of intersection
    … …  
    767766      if (v.z() > 0)
    768767      {
    769         tolORMin  = fRmin1 - 0.5*kRadTolerance*secRMin ;
    770         tolIRMin  = fRmin1 + 0.5*kRadTolerance*secRMin ;
    771         tolIRMax  = fRmax1 - 0.5*kRadTolerance*secRMin ;
    772         tolORMax2 = (fRmax1 + 0.5*kRadTolerance*secRMax)*
    773                     (fRmax1 + 0.5*kRadTolerance*secRMax) ;
     768        tolORMin  = fRmin1 - halfRadTolerance*secRMin ;
     769        tolIRMin  = fRmin1 + halfRadTolerance*secRMin ;
     770        tolIRMax  = fRmax1 - halfRadTolerance*secRMin ;
     771        tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
     772                    (fRmax1 + halfRadTolerance*secRMax) ;
    774773      }
    775774      else
    776775      {
    777         tolORMin  = fRmin2 - 0.5*kRadTolerance*secRMin ;
    778         tolIRMin  = fRmin2 + 0.5*kRadTolerance*secRMin ;
    779         tolIRMax  = fRmax2 - 0.5*kRadTolerance*secRMin ;
    780         tolORMax2 = (fRmax2 + 0.5*kRadTolerance*secRMax)*
    781                     (fRmax2 + 0.5*kRadTolerance*secRMax) ;
     776        tolORMin  = fRmin2 - halfRadTolerance*secRMin ;
     777        tolIRMin  = fRmin2 + halfRadTolerance*secRMin ;
     778        tolIRMax  = fRmax2 - halfRadTolerance*secRMin ;
     779        tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
     780                    (fRmax2 + halfRadTolerance*secRMax) ;
    782781      }
    783782      if ( tolORMin > 0 )
    … …  
    791790        tolIRMin2 = 0.0 ;
    792791      }
    793       if ( tolIRMax > 0 )   tolIRMax2 = tolIRMax*tolIRMax ;      
    794       else                  tolIRMax2 = 0.0 ;
     792      if ( tolIRMax > 0 )  { tolIRMax2 = tolIRMax*tolIRMax; }     
     793      else                 { tolIRMax2 = 0.0; }
    795794     
    796       if (tolIRMin2 <= rhoi2 && rhoi2 <= tolIRMax2)
    797       {
    798   if ( seg && rhoi2 )
    799   {
    800     // Psi = angle made with central (average) phi of shape
    801 
    802     cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    803 
    804     if (cosPsi >= cosHDPhiIT) return s ;
    805   }
    806   else return s ;
    807       }
    808       /*
    809       else if (tolORMin2 <= rhoi2 && rhoi2 <= tolORMax2)
    810       {
    811   if ( seg && rhoi2 )
    812   {
    813     // Psi = angle made with central (average) phi of shape
    814 
    815     cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    816 
    817     if (cosPsi >= cosHDPhiIT) return s ;
    818   }
    819   else return s ;
    820       }
    821       */
     795      if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
     796      {
     797        if ( !fPhiFullCone && rhoi2 )
     798        {
     799          // Psi = angle made with central (average) phi of shape
     800
     801          cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     802
     803          if (cosPsi >= cosHDPhiIT)  { return s; }
     804        }
     805        else
     806        {
     807          return s;
     808        }
     809      }
    822810    }
    823811    else  // On/outside extent, and heading away  -> cannot intersect
    … …  
    861849  if (std::fabs(nt1) > kRadTolerance)  // Equation quadratic => 2 roots
    862850  {
    863       b = nt2/nt1 ;
    864       c = nt3/nt1 ;
    865       d = b*b-c   ;
    866     if ( nt3 > rout*kRadTolerance*secRMax || rout < 0 )
     851    b = nt2/nt1;
     852    c = nt3/nt1;
     853    d = b*b-c  ;
     854    if ( (nt3 > rout*kRadTolerance*secRMax) || (rout < 0) )
    867855    {
    868856      // If outside real cone (should be rho-rout>kRadTolerance*0.5
    869857      // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
    870858
    871 
    872859      if (d >= 0)
    873860      {
    874861         
    875         if (rout < 0 && nt3 <= 0 )
     862        if ((rout < 0) && (nt3 <= 0))
    876863        {
    877864          // Inside `shadow cone' with -ve radius
    … …  
    882869        else
    883870        {
    884           if ( b <= 0 && c >= 0 ) // both >=0, try smaller root
     871          if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
    885872          {
    886873            s = -b - std::sqrt(d) ;
    … …  
    906893            // Z ok. Check phi intersection if reqd
    907894
    908             if ( ! seg )  return s ;
     895            if ( fPhiFullCone )  { return s; }
    909896            else
    910897            {
    … …  
    914901              cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    915902
    916               if ( cosPsi >= cosHDPhiIT ) return s ;
     903              if ( cosPsi >= cosHDPhiIT )  { return s; }
    917904            }
    918905          }
    … …  
    925912      // check not inside, and heading through G4Cons (-> 0 to in)
    926913
    927       if ( t3  > (rin + kRadTolerance*0.5*secRMin)*
    928                  (rin + kRadTolerance*0.5*secRMin) &&
    929            nt2 < 0                                 &&
    930            d >= 0                                  &&
    931            //  nt2 < -kCarTolerance*secRMax/2/fDz                  &&
    932            // t2 < std::sqrt(t3)*v.z()*tanRMax                          &&
    933            // d > kCarTolerance*secRMax*(rout-b*tanRMax*v.z())/nt1 &&
    934            std::fabs(p.z()) <= tolIDz )
     914      if ( ( t3  > (rin + halfRadTolerance*secRMin)*
     915                   (rin + halfRadTolerance*secRMin) )
     916        && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
    935917      {
    936918        // Inside cones, delta r -ve, inside z extent
    937919
    938         if (seg)
     920        if ( !fPhiFullCone )
    939921        {
    940922          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
    941923
    942           if (cosPsi >= cosHDPhiIT) return 0.0 ;
    943         }
    944         else  return 0.0 ;
     924          if (cosPsi >= cosHDPhiIT)  { return 0.0; }
     925        }
     926        else  { return 0.0; }
    945927      }
    946928    }
    … …  
    952934      s = -0.5*nt3/nt2 ;
    953935
    954       if ( s < 0 )   return kInfinity ;   // travel away
     936      if ( s < 0 )  { return kInfinity; }   // travel away
    955937      else  // s >= 0,  If 'forwards'. Check z intersection
    956938      {
    957939        zi = p.z() + s*v.z() ;
    958940
    959         if (std::fabs(zi) <= tolODz && nt2 < 0)
     941        if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
    960942        {
    961943          // Z ok. Check phi intersection if reqd
    962944
    963           if ( ! seg ) return s ;
     945          if ( fPhiFullCone )  { return s; }
    964946          else
    965947          {
    … …  
    969951            cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    970952
    971             if (cosPsi >= cosHDPhiIT) return s ;
     953            if (cosPsi >= cosHDPhiIT)  { return s; }
    972954          }
    973955        }
    … …  
    1015997            if ( std::fabs(zi) <= tolODz )
    1016998            {
    1017               if ( seg )
     999              if ( !fPhiFullCone )
    10181000              {
    10191001                xi     = p.x() + s*v.x() ;
    … …  
    10221004                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    10231005
    1024                 if (cosPsi >= cosHDPhiIT) snxt = s ;
     1006                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
    10251007              }
    1026               else return s ;
     1008              else  { return s; }
    10271009            }
    10281010          }
    … …  
    10461028          ri = rMinAv + zi*tanRMin ;
    10471029
    1048           if ( ri >= 0 )
     1030          if ( ri > 0 )
    10491031          {
    1050             if ( s >= 0 && std::fabs(zi) <= tolODz )  // s > 0
     1032            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s > 0
    10511033            {
    1052               if ( seg )
     1034              if ( !fPhiFullCone )
    10531035              {
    10541036                xi     = p.x() + s*v.x() ;
    … …  
    10561038                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    10571039
    1058                 if (cosPsi >= cosHDPhiOT) snxt = s ;
     1040                if (cosPsi >= cosHDPhiOT)  { snxt = s; }
    10591041              }
    1060               else  return s ;
     1042              else  { return s; }
    10611043            }
    10621044          }
    … …  
    10671049            ri = rMinAv + zi*tanRMin ;
    10681050
    1069             if ( s >= 0 && ri >= 0 && std::fabs(zi) <= tolODz ) // s>0
     1051            if ( (s >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // s>0
    10701052            {
    1071               if ( seg )
     1053              if ( !fPhiFullCone )
    10721054              {
    10731055                xi     = p.x() + s*v.x() ;
    … …  
    10751057                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    10761058
    1077                 if (cosPsi >= cosHDPhiIT) snxt = s ;
     1059                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
    10781060              }
    1079               else  return s ;
     1061              else  { return s; }
    10801062            }
    10811063          }
    … …  
    10951077            // Inside inner real cone, heading outwards, inside z range
    10961078
    1097             if ( seg )
     1079            if ( !fPhiFullCone )
    10981080            {
    10991081              cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
    11001082
    1101               if (cosPsi >= cosHDPhiIT)  return 0.0 ;
     1083              if (cosPsi >= cosHDPhiIT)  { return 0.0; }
    11021084            }
    1103             else  return 0.0 ;
     1085            else  { return 0.0; }
    11041086          }
    11051087          else
    … …  
    11231105                zi = p.z() + s*v.z() ;
    11241106
    1125                 if ( s >= 0 && std::fabs(zi) <= tolODz )  // s>0
     1107                if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
    11261108                {
    1127                   if ( seg )
     1109                  if ( !fPhiFullCone )
    11281110                  {
    11291111                    xi     = p.x() + s*v.x() ;
    … …  
    11321114                    cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    11331115
    1134                     if ( cosPsi >= cosHDPhiIT )  snxt = s ;
     1116                    if ( cosPsi >= cosHDPhiIT )  { snxt = s; }
    11351117                  }
    1136                   else return s ;
     1118                  else  { return s; }
    11371119                }
    11381120              }
    1139               else  return kInfinity ;
     1121              else  { return kInfinity; }
    11401122            }
    11411123          }
    … …  
    11521134            zi = p.z() + s*v.z() ;
    11531135
    1154             if ( s >= 0 && std::fabs(zi) <= tolODz )  // s>0
     1136            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
    11551137            {
    1156               if ( seg )
     1138              if ( !fPhiFullCone )
    11571139              {
    11581140                xi     = p.x() + s*v.x();
    … …  
    11611143                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
    11621144
    1163                 if (cosPsi >= cosHDPhiIT) snxt = s ;
     1145                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
    11641146              }
    1165               else  return s ;
     1147              else  { return s; }
    11661148            }
    11671149          }
    … …  
    11801162  //         -> Should use some form of loop Construct
    11811163
    1182   if ( seg )
    1183   {
    1184     // First phi surface (`S'tarting phi)
    1185 
    1186     sinSPhi = std::sin(fSPhi) ;
    1187     cosSPhi = std::cos(fSPhi) ;
     1164  if ( !fPhiFullCone )
     1165  {
     1166    // First phi surface (starting phi)
     1167
    11881168    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
    11891169                   
    … …  
    11921172      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
    11931173
    1194       if (Dist < kCarTolerance*0.5)
     1174      if (Dist < halfCarTolerance)
    11951175      {
    11961176        s = Dist/Comp ;
    … …  
    11981178        if ( s < snxt )
    11991179        {
    1200           if ( s < 0 ) s = 0.0 ;
     1180          if ( s < 0 )  { s = 0.0; }
    12011181
    12021182          zi = p.z() + s*v.z() ;
    … …  
    12101190            tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
    12111191
    1212             if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2 )
     1192            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
    12131193            {
    12141194              // z and r intersections good - check intersecting with
    12151195              // correct half-plane
    12161196
    1217               if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) snxt = s ;
    1218             }   
     1197              if ((yi*cosCPhi - xi*sinCPhi) <= 0 )  { snxt = s; }
     1198            }
    12191199          }
    12201200        }
    12211201      }
    1222     }   
    1223     // Second phi surface (`E'nding phi)
    1224 
    1225     ePhi    = fSPhi + fDPhi ;
    1226     sinEPhi = std::sin(ePhi) ;
    1227     cosEPhi = std::cos(ePhi) ;
     1202    }
     1203
     1204    // Second phi surface (Ending phi)
     1205
    12281206    Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
    12291207       
    … …  
    12311209    {
    12321210      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
    1233       if (Dist < kCarTolerance*0.5)
     1211      if (Dist < halfCarTolerance)
    12341212      {
    12351213        s = Dist/Comp ;
    … …  
    12371215        if ( s < snxt )
    12381216        {
    1239           if ( s < 0 )  s = 0.0 ;
     1217          if ( s < 0 )  { s = 0.0; }
    12401218
    12411219          zi = p.z() + s*v.z() ;
    … …  
    12491227            tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
    12501228
    1251             if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2 )
     1229            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
    12521230            {
    12531231              // z and r intersections good - check intersecting with
    12541232              // correct half-plane
    12551233
    1256               if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  snxt = s ;
    1257             }   
     1234              if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  { snxt = s; }
     1235            }
    12581236          }
    12591237        }
    … …  
    12611239    }
    12621240  }
    1263   if (snxt < kCarTolerance*0.5) snxt = 0.;
    1264 
    1265 #ifdef consdebug
    1266   G4cout.precision(24);
    1267   G4cout<<"G4Cons::DistanceToIn(p,v) "<<G4endl;
    1268   G4cout<<"position = "<<p<<G4endl;
    1269   G4cout<<"direction = "<<v<<G4endl;
    1270   G4cout<<"distance = "<<snxt<<G4endl;
    1271 #endif
     1241  if (snxt < halfCarTolerance)  { snxt = 0.; }
    12721242
    12731243  return snxt ;
    … …  
    12831253G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
    12841254{
    1285   G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
     1255  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
    12861256  G4double tanRMin, secRMin, pRMin ;
    12871257  G4double tanRMax, secRMax, pRMax ;
    1288   G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi ;
    1289   G4double cosPsi ;
    12901258
    12911259  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
    … …  
    13041272    safeR2  = (rho - pRMax)/secRMax ;
    13051273
    1306     if ( safeR1 > safeR2) safe = safeR1 ;
    1307     else                  safe = safeR2 ;
     1274    if ( safeR1 > safeR2) { safe = safeR1; }
     1275    else                  { safe = safeR2; }
    13081276  }
    13091277  else
    … …  
    13141282    safe    = (rho - pRMax)/secRMax ;
    13151283  }
    1316   if ( safeZ > safe ) safe = safeZ ;
    1317 
    1318   if ( fDPhi < twopi && rho )
    1319   {
    1320     phiC    = fSPhi + fDPhi*0.5 ;
    1321     cosPhiC = std::cos(phiC) ;
    1322     sinPhiC = std::sin(phiC) ;
    1323 
     1284  if ( safeZ > safe )  { safe = safeZ; }
     1285
     1286  if ( !fPhiFullCone && rho )
     1287  {
    13241288    // Psi=angle from central phi to point
    13251289
    1326     cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
     1290    cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
    13271291
    13281292    if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
    13291293    {
    1330       if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0.0 )
    1331       {
    1332         safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
     1294      if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
     1295      {
     1296        safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
    13331297      }
    13341298      else
    13351299      {
    1336         ePhi    = fSPhi + fDPhi ;
    1337         safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
    1338       }
    1339       if ( safePhi > safe ) safe = safePhi ;
    1340     }
    1341   }
    1342   if ( safe < 0.0 ) safe = 0.0 ;
     1300        safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
     1301      }
     1302      if ( safePhi > safe )  { safe = safePhi; }
     1303    }
     1304  }
     1305  if ( safe < 0.0 )  { safe = 0.0; }
    13431306
    13441307  return safe ;
    … …  
    13471310///////////////////////////////////////////////////////////////
    13481311//
    1349 // Calculate distance to surface of shape from `inside', allowing for tolerance
     1312// Calculate distance to surface of shape from 'inside', allowing for tolerance
    13501313// - Only Calc rmax intersection if no valid rmin intersection
    13511314
    13521315G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
    1353               const G4ThreeVector& v,
    1354               const G4bool calcNorm,
    1355                     G4bool *validNorm,
    1356                     G4ThreeVector *n) const
     1316                                const G4ThreeVector& v,
     1317                                const G4bool calcNorm,
     1318                                      G4bool *validNorm,
     1319                                      G4ThreeVector *n) const
    13571320{
    13581321  ESide side = kNull, sider = kNull, sidephi = kNull;
    13591322
     1323  static const G4double halfCarTolerance=kCarTolerance*0.5;
     1324  static const G4double halfRadTolerance=kRadTolerance*0.5;
     1325  static const G4double halfAngTolerance=kAngTolerance*0.5;
     1326
    13601327  G4double snxt,sr,sphi,pdist ;
    13611328
    … …  
    13731340  // Vars for phi intersection:
    13741341
    1375   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
    1376   G4double cPhi, sinCPhi, cosCPhi ;
    13771342  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
    13781343  G4double zi, ri, deltaRoi2 ;
    … …  
    13841349    pdist = fDz - p.z() ;
    13851350
    1386     if (pdist > kCarTolerance*0.5)
     1351    if (pdist > halfCarTolerance)
    13871352    {
    13881353      snxt = pdist/v.z() ;
    … …  
    13961361        *validNorm = true ;
    13971362      }
    1398       return snxt = 0.0 ;
     1363      return  snxt = 0.0;
    13991364    }
    14001365  }
    … …  
    14031368    pdist = fDz + p.z() ;
    14041369
    1405     if ( pdist > kCarTolerance*0.5)
     1370    if ( pdist > halfCarTolerance)
    14061371    {
    14071372      snxt = -pdist/v.z() ;
    … …  
    14651430                - fRmax1*(fRmax1 + kRadTolerance*secRMax);
    14661431  }
    1467   else deltaRoi2 = 1.0 ;
    1468 
    1469   if ( nt1 && deltaRoi2 > 0.0 ) 
     1432  else
     1433  {
     1434    deltaRoi2 = 1.0;
     1435  }
     1436
     1437  if ( nt1 && (deltaRoi2 > 0.0) ) 
    14701438  {
    14711439    // Equation quadratic => 2 roots : second root must be leaving
    … …  
    14781446    {
    14791447      // Check if on outer cone & heading outwards
    1480       // NOTE: Should use rho-rout>-kRadtolerance*0.5
     1448      // NOTE: Should use rho-rout>-kRadTolerance*0.5
    14811449       
    1482       if (nt3 > -kRadTolerance*0.5 && nt2 >= 0 )
     1450      if (nt3 > -halfRadTolerance && nt2 >= 0 )
    14831451      {
    14841452        if (calcNorm)
    … …  
    14861454          risec      = std::sqrt(t3)*secRMax ;
    14871455          *validNorm = true ;
    1488           *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
     1456          *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
    14891457        }
    14901458        return snxt=0 ;
    … …  
    14971465        ri    = tanRMax*zi + rMaxAv ;
    14981466         
    1499         if ( (ri >= 0) && (-kRadTolerance*0.5 <= sr) &&
    1500                           ( sr <= kRadTolerance*0.5)     )
     1467        if ((ri >= 0) && (-halfRadTolerance <= sr) && (sr <= halfRadTolerance))
    15011468        {
    15021469          // An intersection within the tolerance
    … …  
    15061473          sidetol = kRMax ;
    15071474        }           
    1508         if ( (ri < 0) || (sr < kRadTolerance*0.5) )
     1475        if ( (ri < 0) || (sr < halfRadTolerance) )
    15091476        {
    15101477          // Safety: if both roots -ve ensure that sr cannot `win'
    … …  
    15151482          ri  = tanRMax*zi + rMaxAv ;
    15161483
    1517           if (ri >= 0 && sr2 > kRadTolerance*0.5)  sr = sr2 ;
     1484          if ((ri >= 0) && (sr2 > halfRadTolerance))
     1485          {
     1486            sr = sr2;
     1487          }
    15181488          else
    15191489          {
    15201490            sr = kInfinity ;
    15211491
    1522             if(    (-kRadTolerance*0.5 <= sr2)
    1523                 && ( sr2 <= kRadTolerance*0.5)  )
     1492            if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
    15241493            {
    15251494              // An intersection within the tolerance.
    … …  
    15401509      if ( calcNorm )
    15411510      {
    1542         risec      = std::sqrt(t3)*secRMax ;
    1543         *validNorm = true ;
    1544         *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
     1511        risec      = std::sqrt(t3)*secRMax;
     1512        *validNorm = true;
     1513        *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
    15451514      }
    15461515      return snxt = 0.0 ;
    15471516    }
    15481517  }
    1549   else if ( nt2 && deltaRoi2 > 0.0 )
     1518  else if ( nt2 && (deltaRoi2 > 0.0) )
    15501519  {
    15511520    // Linear case (only one intersection) => point outside outer cone
    … …  
    15531522    if ( calcNorm )
    15541523    {
    1555       risec      = std::sqrt(t3)*secRMax ;
    1556       *validNorm = true ;
    1557       *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
     1524      risec      = std::sqrt(t3)*secRMax;
     1525      *validNorm = true;
     1526      *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
    15581527    }
    15591528    return snxt = 0.0 ;
    … …  
    15691538  // Check possible intersection within tolerance
    15701539
    1571   if ( slentol <= kCarTolerance*0.5 )
     1540  if ( slentol <= halfCarTolerance )
    15721541  {
    15731542    // An intersection within the tolerance was found. 
    … …  
    15801549    // Calculate a normal vector,  as below
    15811550
    1582     xi    = p.x() + slentol*v.x() ;
    1583     yi    = p.y() + slentol*v.y() ;
    1584     risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
    1585     G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
     1551    xi    = p.x() + slentol*v.x();
     1552    yi    = p.y() + slentol*v.y();
     1553    risec = std::sqrt(xi*xi + yi*yi)*secRMax;
     1554    G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
    15861555
    15871556    if ( Normal.dot(v) > 0 )    // We will leave the Cone immediatelly
    … …  
    15951564    }
    15961565    else // On the surface, but not heading out so we ignore this intersection
    1597     {    //                                    (as it is within tolerance). 
     1566    {    //                                        (as it is within tolerance).
    15981567      slentol = kInfinity ;
    15991568    }
    … …  
    16211590      d = b*b - c ;
    16221591
    1623       if (d >= 0.0 )
     1592      if ( d >= 0.0 )
    16241593      {
    16251594        // NOTE: should be rho-rin<kRadTolerance*0.5,
    … …  
    16301599          if ( nt2 < 0.0 )
    16311600          {
    1632             if (calcNorm)  *validNorm = false ;
    1633             return          snxt      = 0.0 ;
     1601            if (calcNorm)  { *validNorm = false; }
     1602            return          snxt      = 0.0;
    16341603          }
    16351604        }
    … …  
    16401609          ri  = tanRMin*zi + rMinAv ;
    16411610
    1642           if( (ri >= 0.0) && (-kRadTolerance*0.5 <= sr2) &&
    1643                              ( sr2 <= kRadTolerance*0.5)      )
     1611          if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
    16441612          {
    16451613            // An intersection within the tolerance
    … …  
    16491617            sidetol = kRMax ;
    16501618          }
    1651           if( (ri<0) || (sr2 < kRadTolerance*0.5) )
     1619          if( (ri<0) || (sr2 < halfRadTolerance) )
    16521620          {
    16531621            sr3 = -b + std::sqrt(d) ;
    … …  
    16561624            //         distancetoout
    16571625
    1658             if  ( sr3 > kCarTolerance*0.5 )
     1626            if  ( sr3 > halfRadTolerance )
    16591627            {
    16601628              if( sr3 < sr )
    … …  
    16701638              }
    16711639            }
    1672             else if ( sr3 > -kCarTolerance*0.5 )
     1640            else if ( sr3 > -halfRadTolerance )
    16731641            {
    16741642              // Intersection in tolerance. Store to check if it's good
    … …  
    16781646            }
    16791647          }
    1680           else if ( sr2 < sr && sr2 > kCarTolerance*0.5 )
     1648          else if ( (sr2 < sr) && (sr2 > halfCarTolerance) )
    16811649          {
    16821650            sr    = sr2 ;
    16831651            sider = kRMin ;
    16841652          }
    1685           else if (sr2 > -kCarTolerance*0.5)
     1653          else if (sr2 > -halfCarTolerance)
    16861654          {
    16871655            // Intersection in tolerance. Store to check if it's good
    … …  
    16901658            sidetol = kRMin ;
    16911659          }   
    1692           if( slentol <= kCarTolerance*0.5  )
     1660          if( slentol <= halfCarTolerance  )
    16931661          {
    16941662            // An intersection within the tolerance was found.
    … …  
    17061674            if( Normal.dot(v) > 0 )
    17071675            {
    1708               // We will leave the Cone immediatelly
     1676              // We will leave the cone immediately
     1677
    17091678              if( calcNorm )
    17101679              {
    … …  
    17311700  // Phi Intersection
    17321701 
    1733   if ( fDPhi < twopi )
    1734   {
    1735     sinSPhi = std::sin(fSPhi) ;
    1736     cosSPhi = std::cos(fSPhi) ;
    1737     ePhi    = fSPhi + fDPhi ;
    1738     sinEPhi = std::sin(ePhi) ;
    1739     cosEPhi = std::cos(ePhi) ;
    1740     cPhi    = fSPhi + fDPhi*0.5 ;
    1741     sinCPhi = std::sin(cPhi) ;
    1742     cosCPhi = std::cos(cPhi) ;
     1702  if ( !fPhiFullCone )
     1703  {
    17431704    // add angle calculation with correction
    1744       // of the  difference in domain of atan2 and Sphi
    1745         vphi = std::atan2(v.y(),v.x()) ;
    1746 
    1747          if ( vphi < fSPhi - kAngTolerance*0.5  )             vphi += twopi ;
    1748          else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) vphi -= twopi;
     1705    // of the difference in domain of atan2 and Sphi
     1706
     1707    vphi = std::atan2(v.y(),v.x()) ;
     1708
     1709    if ( vphi < fSPhi - halfAngTolerance  )              { vphi += twopi; }
     1710    else if ( vphi > fSPhi + fDPhi + halfAngTolerance )  { vphi -= twopi; }
     1711
    17491712    if ( p.x() || p.y() )   // Check if on z axis (rho not needed later)
    17501713    {
    … …  
    17611724      sidephi = kNull ;
    17621725
    1763       if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
    1764                               && (pDistE <= 0.5*kCarTolerance) ) )
    1765          || ( (fDPhi >  pi) && !((pDistS >  0.5*kCarTolerance)
    1766                               && (pDistE >  0.5*kCarTolerance) ) )  )
    1767         {
    1768           // Inside both phi *full* planes
    1769           if ( compS < 0 )
     1726      if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
     1727                            && (pDistE <= halfCarTolerance) ) )
     1728         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
     1729                              && (pDistE >  halfCarTolerance) ) )  )
     1730      {
     1731        // Inside both phi *full* planes
     1732        if ( compS < 0 )
     1733        {
     1734          sphi = pDistS/compS ;
     1735          if (sphi >= -halfCarTolerance)
    17701736          {
    1771             sphi = pDistS/compS ;
    1772             if (sphi >= -0.5*kCarTolerance)
     1737            xi = p.x() + sphi*v.x() ;
     1738            yi = p.y() + sphi*v.y() ;
     1739
     1740            // Check intersecting with correct half-plane
     1741            // (if not -> no intersect)
     1742            //
     1743            if ( (std::abs(xi)<=kCarTolerance)
     1744              && (std::abs(yi)<=kCarTolerance) )
    17731745            {
    1774               xi = p.x() + sphi*v.x() ;
    1775               yi = p.y() + sphi*v.y() ;
    1776 
    1777               // Check intersecting with correct half-plane
    1778               // (if not -> no intersect)
    1779               //
    1780                if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){
    1781                  sidephi= kSPhi;
    1782                 if(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi))
    1783                     { sphi = kInfinity; }
    1784                                        
     1746              sidephi= kSPhi;
     1747              if ( ( fSPhi-halfAngTolerance <= vphi )
     1748                && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
     1749              {
     1750                sphi = kInfinity;
    17851751              }
    1786               else
    1787               if ((yi*cosCPhi-xi*sinCPhi)>=0)
    1788               {
    1789                 sphi = kInfinity ;
    1790               }
    1791               else
    1792               {
    1793                 sidephi = kSPhi ;
    1794                 if ( pDistS > -kCarTolerance*0.5 )
    1795                 {
    1796                   sphi = 0.0 ; // Leave by sphi immediately
    1797                 }   
    1798               }       
     1752            }
     1753            else
     1754            if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
     1755            {
     1756              sphi = kInfinity ;
    17991757            }
    18001758            else
    18011759            {
    1802               sphi = kInfinity ;
    1803             }
     1760              sidephi = kSPhi ;
     1761              if ( pDistS > -halfCarTolerance )
     1762              {
     1763                sphi = 0.0 ; // Leave by sphi immediately
     1764              }   
     1765            }       
    18041766          }
    18051767          else
    … …  
    18071769            sphi = kInfinity ;
    18081770          }
    1809 
    1810           if ( compE < 0 )
     1771        }
     1772        else
     1773        {
     1774          sphi = kInfinity ;
     1775        }
     1776
     1777        if ( compE < 0 )
     1778        {
     1779          sphi2 = pDistE/compE ;
     1780
     1781          // Only check further if < starting phi intersection
     1782          //
     1783          if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
    18111784          {
    1812             sphi2 = pDistE/compE ;
    1813 
    1814             // Only check further if < starting phi intersection
    1815             //
    1816             if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
     1785            xi = p.x() + sphi2*v.x() ;
     1786            yi = p.y() + sphi2*v.y() ;
     1787
     1788            // Check intersecting with correct half-plane
     1789
     1790            if ( (std::abs(xi)<=kCarTolerance)
     1791              && (std::abs(yi)<=kCarTolerance) )
    18171792            {
    1818               xi = p.x() + sphi2*v.x() ;
    1819               yi = p.y() + sphi2*v.y() ;
    1820 
    1821               // Check intersecting with correct half-plane
    1822                if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){
    1823                // Leaving via ending phi
    1824                 if(!(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi))){
    1825                   sidephi = kEPhi ;
    1826                   if ( pDistE <= -kCarTolerance*0.5 ) sphi = sphi2 ;
    1827                   else                                sphi = 0.0 ;
    1828                   }
    1829                 }
    1830              else // Check intersecting with correct half-plane
    1831               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
     1793              // Leaving via ending phi
     1794
     1795              if(!( (fSPhi-halfAngTolerance <= vphi)
     1796                 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
    18321797              {
    1833                 // Leaving via ending phi
    1834 
    18351798                sidephi = kEPhi ;
    1836                 if ( pDistE <= -kCarTolerance*0.5 ) sphi = sphi2 ;
    1837                 else                                sphi = 0.0 ;
     1799                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
     1800                else                                { sphi = 0.0; }
    18381801              }
    18391802            }
     1803            else // Check intersecting with correct half-plane
     1804            if ( yi*cosCPhi-xi*sinCPhi >= 0 )
     1805            {
     1806              // Leaving via ending phi
     1807
     1808              sidephi = kEPhi ;
     1809              if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
     1810              else                                { sphi = 0.0; }
     1811            }
    18401812          }
    18411813        }
    1842         else
    1843         {
    1844           sphi = kInfinity ;
    1845         }
    1846 
    1847 
     1814      }
     1815      else
     1816      {
     1817        sphi = kInfinity ;
     1818      }
    18481819    }
    18491820    else
    … …  
    18521823      // within phi of shape, Step limited by rmax, else Step =0
    18531824
    1854       // vphi = std::atan2(v.y(),v.x()) ;
    1855 
    1856       // if ( fSPhi < vphi && vphi < fSPhi + fDPhi ) sphi = kInfinity ;
    1857 
    1858        if ( ((fSPhi-0.5*kAngTolerance) <= vphi) && (vphi <=( fSPhi + fDPhi)+0.5*kAngTolerance) )
    1859          {
    1860           sphi = kInfinity ;
    1861         }
     1825      if ( (fSPhi-halfAngTolerance <= vphi)
     1826        && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
     1827      {
     1828        sphi = kInfinity ;
     1829      }
    18621830      else
    18631831      {
    … …  
    18681836    if ( sphi < snxt )  // Order intersecttions
    18691837    {
    1870         snxt=sphi ;
    1871         side=sidephi ;
     1838      snxt=sphi ;
     1839      side=sidephi ;
    18721840    }
    18731841  }
    … …  
    18801848  {
    18811849    switch(side)
    1882     {
    1883       case kRMax:
    1884           // Note: returned vector not normalised
    1885           // (divide by frmax for unit vector)
     1850    {                     // Note: returned vector not normalised
     1851      case kRMax:         // (divide by frmax for unit vector)
    18861852        xi         = p.x() + snxt*v.x() ;
    18871853        yi         = p.y() + snxt*v.y() ;
    … …  
    18911857        break ;
    18921858      case kRMin:
    1893         *validNorm=false ;  // Rmin is inconvex
     1859        *validNorm = false ;  // Rmin is inconvex
    18941860        break ;
    18951861      case kSPhi:
    18961862        if ( fDPhi <= pi )
    18971863        {
    1898           *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
     1864          *n         = G4ThreeVector(sinSPhi, -cosSPhi, 0);
    18991865          *validNorm = true ;
    19001866        }
    1901         else   *validNorm = false ;
     1867        else
     1868        {
     1869          *validNorm = false ;
     1870        }
    19021871        break ;
    19031872      case kEPhi:
    19041873        if ( fDPhi <= pi )
    19051874        {
    1906           *n         = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
     1875          *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
    19071876          *validNorm = true ;
    19081877        }
    1909         else  *validNorm = false ;
     1878        else
     1879        {
     1880          *validNorm = false ;
     1881        }
    19101882        break ;
    19111883      case kPZ:
    … …  
    19251897        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
    19261898        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
    1927         G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm"
    1928                << G4endl << G4endl ;
     1899        G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
     1900               << " mm" << G4endl << G4endl ;
    19291901        if( p.x() != 0. || p.x() != 0.)
    19301902        {
    1931            G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree << " degree"
    1932                   << G4endl << G4endl ;
     1903           G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
     1904                  << " degree" << G4endl << G4endl ;
    19331905        }
    19341906        G4cout << "Direction:" << G4endl << G4endl ;
    … …  
    19431915    }
    19441916  }
    1945   if (snxt < kCarTolerance*0.5) snxt = 0.;
    1946 #ifdef consdebug
    1947   G4cout.precision(24);
    1948   G4cout<<"G4Cons::DistanceToOut(p,v,...) "<<G4endl;
    1949   G4cout<<"position = "<<p<<G4endl;
    1950   G4cout<<"direction = "<<v<<G4endl;
    1951   G4cout<<"distance = "<<snxt<<G4endl;
    1952 #endif
     1917  if (snxt < halfCarTolerance)  { snxt = 0.; }
     1918
    19531919  return snxt ;
    19541920}
    … …  
    19601926G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
    19611927{
    1962   G4double safe=0.0,rho,safeR1,safeR2,safeZ ;
    1963   G4double tanRMin,secRMin,pRMin ;
    1964   G4double tanRMax,secRMax,pRMax ;
    1965   G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi ;
     1928  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
     1929  G4double tanRMin, secRMin, pRMin;
     1930  G4double tanRMax, secRMax, pRMax;
    19661931
    19671932#ifdef G4CSGDEBUG
    … …  
    19751940    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
    19761941    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
    1977     G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm"
    1978            << G4endl << G4endl ;
    1979     if( p.x() != 0. || p.x() != 0.)
    1980     {
    1981       G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree << " degree"
    1982              << G4endl << G4endl ;
    1983     }
    1984     G4Exception("G4Cons::DistanceToOut(p)", "Notification", JustWarning,
    1985                 "Point p is outside !?" );
     1942    G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
     1943           << " mm" << G4endl << G4endl ;
     1944    if( (p.x() != 0.) || (p.x() != 0.) )
     1945    {
     1946      G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
     1947             << " degree" << G4endl << G4endl ;
     1948    }
     1949    G4Exception("G4Cons::DistanceToOut(p)", "Notification",
     1950                JustWarning, "Point p is outside !?" );
    19861951  }
    19871952#endif
    … …  
    19971962    safeR1  = (rho - pRMin)/secRMin ;
    19981963  }
    1999   else safeR1 = kInfinity ;
     1964  else
     1965  {
     1966    safeR1 = kInfinity ;
     1967  }
    20001968
    20011969  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
    … …  
    20041972  safeR2  = (pRMax - rho)/secRMax ;
    20051973
    2006   if (safeR1 < safeR2) safe = safeR1 ;
    2007   else                 safe = safeR2 ;
    2008   if (safeZ < safe)    safe = safeZ  ;
     1974  if (safeR1 < safeR2)  { safe = safeR1; }
     1975  else                  { safe = safeR2; }
     1976  if (safeZ < safe)     { safe = safeZ ; }
    20091977
    20101978  // Check if phi divided, Calc distances closest phi plane
    20111979
    2012   if (fDPhi < twopi)
     1980  if (!fPhiFullCone)
    20131981  {
    20141982    // Above/below central phi of G4Cons?
    20151983
    2016     phiC    = fSPhi + fDPhi*0.5 ;
    2017     cosPhiC = std::cos(phiC) ;
    2018     sinPhiC = std::sin(phiC) ;
    2019 
    2020     if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
    2021     {
    2022       safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
     1984    if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
     1985    {
     1986      safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
    20231987    }
    20241988    else
    20251989    {
    2026       ePhi    = fSPhi + fDPhi ;
    2027       safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
    2028     }
    2029     if (safePhi < safe) safe = safePhi ;
    2030   }
    2031   if ( safe < 0 ) safe = 0 ;
    2032   return safe ; 
     1990      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
     1991    }
     1992    if (safePhi < safe)  { safe = safePhi; }
     1993  }
     1994  if ( safe < 0 )  { safe = 0; }
     1995
     1996  return safe ;
    20331997}
    20341998
    … …  
    20682032  meshAngle = fDPhi/(noCrossSections - 1) ;
    20692033
    2070   // G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1  ;
    2071 
    20722034  meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
    20732035  meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
    … …  
    20762038  // on the x axis. Will give better extent calculations when not rotated.
    20772039
    2078   if (fDPhi == twopi && fSPhi == 0.0 )
     2040  if ( fPhiFullCone && (fSPhi == 0.0) )
    20792041  {
    20802042    sAngle = -meshAngle*0.5 ;
    … …  
    21022064      rMaxY2 = meshRMax2*sinCrossAngle ;
    21032065       
    2104       // G4double RMin = (fRmin2 <= fRmin1) ? fRmin2 : fRmin1  ;
    2105 
    21062066      rMinX1 = fRmin1*cosCrossAngle ;
    21072067      rMinY1 = fRmin1*sinCrossAngle ;
    … …  
    21272087                "Error in allocation of vertices. Out of memory !");
    21282088  }
     2089
    21292090  return vertices ;
    21302091}
    … …  
    21922153  rRand2 = RandFlat::shoot(fRmin2,fRmax2);
    21932154 
    2194   if(fSPhi == 0. && fDPhi == twopi){ Afive = 0.; }
     2155  if ( (fSPhi == 0.) && fPhiFullCone )  { Afive = 0.; }
    21952156  chose  = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
    21962157 
    … …  
    22252186  else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
    22262187  {
    2227     return G4ThreeVector (rRand1*cosu,rRand1*sinu,-1*fDz);
     2188    return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
    22282189  }
    22292190  else if( (chose >= Aone + Atwo + Athree)
Note: See TracChangeset for help on using the changeset viewer.