Ignore:
Timestamp:
Jan 8, 2010, 11:56:51 AM (14 years ago)
Author:
garnier
Message:

update geant4.9.3 tag

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

Legend:

Unmodified
Added
Removed

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

    r921 r1228  
    1 $Id: History,v 1.109 2008/11/21 09:50:20 gcosmo Exp $
     1$Id: History,v 1.119 2009/11/30 10:20:53 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20Nov 30, 2009  V.Grichine                 geom-csg-V09-02-08
     21- G4Orb: moved debug warning in DistanceToIn(p,v) within G4CSGDEBUG flag.
     22
     23Nov 26, 2009  T.Nikitina                 geom-csg-V09-02-07
     24- G4Torus: fix in SolveNumericJT() in order to take in account the difference
     25  in the value of theta for different intervals, [0:pi] or [-pi:0], and for
     26  SPhi in [0:twopi] or [-twopi:0]. Addresses problem report #1086.
     27
     28Nov 12, 2009  T.Nikitina                 geom-csg-V09-02-06
     29- G4Cons: fix to DistanceToIn(p,v), added a check on the direction in case
     30  of point on surface. Fixes a problem of stuck tracks observed in CMS, due
     31  to wrong reply from the solid for points on the inner radius surface base
     32  with direction along the imaginary extension of the cone.
     33- Added test case to internal test testG4Cons2.
     34
     35Aug 07, 2009  T.Nikitina                 geom-csg-V09-02-05
     36- G4Sphere: fix for the calculation of the normal in DistanceToOut()
     37  to avoid cases of division by zero in specific configurations.
     38  Addresses problem report #977.
     39
     40Jun 30, 2009  T.Nikitina                 geom-csg-V09-02-04
     41- Introduced to DistanceToIn(p,v) splitting of the distance for point very far
     42  from intersection area and big difference between solid dimensions and
     43  distance to it; resolves issue observed on 64 bits problem. Affected solids:
     44  G4Tubs, G4Cons, G4Sphere, G4Orb. Addresses problem report #1022.
     45- Fixed rare cases of NaN on G4Sphere.
     46
     47Jun 09, 2009  G.Cosmo                    geom-csg-V09-02-03
     48- Fixed case of non-initialised start-Phi value in G4Tubs and G4Cons,
     49  problem introduced following modifications in tag "geom-csg-V09-02-00".
     50
     51May 14, 2009  T.Nikitina                 geom-csg-V09-02-02
     52- G4Sphere: initialise flags for full sphere in constructor.
     53
     54May 14, 2009  T.Nikitina                 geom-csg-V09-02-01
     55- G4Sphere: correction di DistanceToOut(p,v, ...) for phi sections for
     56  rays passing though zero.
     57- More minor refiniments to G4Sphere following review.
     58
     59May 08, 2009  G.Cosmo                    geom-csg-V09-02-00
     60- G4Sphere: implemented speed improvements and corrections from joint
     61  code review of G4Sphere class. Cached computation for half-tolerances
     62  and use of Boolean flag for identifying if full-sphere, shell or section.
     63  Implemented caching of trigonometric values, now directly computed inside
     64  modifiers for Phi and Theta angles as required for parametrised cases.
     65  Rationalised usage of relative radial tolerances.
     66- G4Tubs, G4Cons: rationalised usage of modifiers for Phi angles and
     67  simplified constructors.
     68- Added new test case for relatively big spheres in unit tests testG4Orb
     69  and testG4Sphere, based on problem report #1046.
    1970
    2071Nov 21, 2008  G.Cosmo                    geom-csg-V09-01-08

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

    r1058 r1228  
    2626//
    2727// $Id: G4Box.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2626//
    2727// $Id: G4Box.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2626//
    2727// $Id: G4CSGSolid.hh,v 1.12 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// 

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Cons.hh,v 1.21 2008/11/06 11:04:00 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Cons.hh,v 1.22 2009/03/31 09:56:24 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    8080                 G4double pDz,
    8181                 G4double pSPhi, G4double pDPhi);
    82          
    83     virtual ~G4Cons() ;
     82      //
     83      // Constructs a cone with the given name and dimensions
     84
     85   ~G4Cons() ;
     86      //
     87      // Destructor
    8488
    8589    // Accessors
    8690
    87     inline G4double    GetInnerRadiusMinusZ() const;
    88     inline G4double    GetOuterRadiusMinusZ() const;
    89     inline G4double    GetInnerRadiusPlusZ()  const;
    90     inline G4double    GetOuterRadiusPlusZ()  const;
    91  
    92     inline G4double    GetZHalfLength()       const;
    93  
    94     inline G4double    GetStartPhiAngle () const;
    95     inline G4double    GetDeltaPhiAngle () const;
     91    inline G4double GetInnerRadiusMinusZ() const;
     92    inline G4double GetOuterRadiusMinusZ() const;
     93    inline G4double GetInnerRadiusPlusZ()  const;
     94    inline G4double GetOuterRadiusPlusZ()  const;
     95    inline G4double GetZHalfLength()       const;
     96    inline G4double GetStartPhiAngle()     const;
     97    inline G4double GetDeltaPhiAngle()     const;
    9698 
    9799    // Modifiers
    98100
    99     inline void    SetInnerRadiusMinusZ( G4double Rmin1 );
    100     inline void    SetOuterRadiusMinusZ( G4double Rmax1 );
    101     inline void    SetInnerRadiusPlusZ ( G4double Rmin2 );
    102     inline void    SetOuterRadiusPlusZ ( G4double Rmax2 );
    103          
    104     inline void    SetZHalfLength      ( G4double newDz );
    105     inline void    SetStartPhiAngle    ( G4double newSPhi);
    106     inline void    SetDeltaPhiAngle    ( G4double newDPhi);
     101    inline void SetInnerRadiusMinusZ (G4double Rmin1 );
     102    inline void SetOuterRadiusMinusZ (G4double Rmax1 );
     103    inline void SetInnerRadiusPlusZ  (G4double Rmin2 );
     104    inline void SetOuterRadiusPlusZ  (G4double Rmax2 );
     105    inline void SetZHalfLength       (G4double newDz );
     106    inline void SetStartPhiAngle     (G4double newSPhi, G4bool trig=true);
     107    inline void SetDeltaPhiAngle     (G4double newDPhi);
    107108
    108109    // Other methods for solid
    109110
    110     inline G4double    GetCubicVolume();
    111     inline G4double    GetSurfaceArea();
    112 
    113     void ComputeDimensions(G4VPVParameterisation* p,
    114                            const G4int n,
    115                            const G4VPhysicalVolume* pRep);
    116 
    117     G4bool CalculateExtent(const EAxis pAxis,
    118                            const G4VoxelLimits& pVoxelLimit,
    119                            const G4AffineTransform& pTransform,
    120                                  G4double& pmin, G4double& pmax) const;         
    121 
    122     EInside Inside(const G4ThreeVector& p) const;
    123 
    124     G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const;
     111    inline G4double GetCubicVolume();
     112    inline G4double GetSurfaceArea();
     113
     114    void ComputeDimensions(       G4VPVParameterisation* p,
     115                            const G4int n,
     116                            const G4VPhysicalVolume* pRep );
     117
     118    G4bool CalculateExtent( const EAxis pAxis,
     119                            const G4VoxelLimits& pVoxelLimit,
     120                            const G4AffineTransform& pTransform,
     121                                  G4double& pmin, G4double& pmax ) const;         
     122
     123    EInside Inside( const G4ThreeVector& p ) const;
     124
     125    G4ThreeVector SurfaceNormal( const G4ThreeVector& p ) const;
    125126
    126127    G4double DistanceToIn (const G4ThreeVector& p,
    … …  
    134135    G4double DistanceToOut(const G4ThreeVector& p) const;
    135136
    136     G4GeometryType  GetEntityType() const;
     137    G4GeometryType GetEntityType() const;
    137138       
    138139    G4ThreeVector GetPointOnSurface() const;
    … …  
    160161    inline G4double    GetRmin2() const;
    161162    inline G4double    GetRmax2() const;
    162  
    163163    inline G4double    GetDz()    const;
    164  
    165164    inline G4double    GetSPhi() const;
    166165    inline G4double    GetDPhi() const;
    167166
    168   protected:
     167  private:
    169168 
    170169    G4ThreeVectorList*
    171170    CreateRotatedVertices(const G4AffineTransform& pTransform) const;
    172171 
    173     G4double fRmin1, fRmin2, fRmax1, fRmax2, fDz, fSPhi, fDPhi;
    174     G4bool fPhiFullCone;
     172    inline void Initialize();
     173      //
     174      // Reset relevant values to zero
     175
     176    inline void CheckSPhiAngle(G4double sPhi);
     177    inline void CheckDPhiAngle(G4double dPhi);
     178    inline void CheckPhiAngles(G4double sPhi, G4double dPhi);
     179      //
     180      // Reset relevant flags and angle values
     181
     182    inline void InitializeTrigonometry();
     183      //
     184      // Recompute relevant trigonometric values and cache them
     185
     186    G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector& p) const;
     187      //
     188      // Algorithm for SurfaceNormal() following the original
     189      // specification for points not on the surface
     190
     191  private:
    175192
    176193    // Used by distanceToOut
    177  
     194    //
    178195    enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
    179196 
    180197    // used by normal
    181  
     198    //
    182199    enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
    183200
    184   private:
    185 
    186     inline void Initialise();
    187       // Reset relevant values to zero
    188 
    189     inline void InitializeTrigonometry();
    190       //
    191       // Recompute relevant trigonometric values and cache them
    192 
    193     G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector& p) const;
    194       //
    195       // Algorithm for SurfaceNormal() following the original
    196       // specification for points not on the surface
    197 
    198   private:
    199 
    200201    G4double kRadTolerance, kAngTolerance;
    201202      //
    202203      // Radial and angular tolerances
     204
     205    G4double fRmin1, fRmin2, fRmax1, fRmax2, fDz, fSPhi, fDPhi;
     206      //
     207      // Radial and angular dimensions
    203208
    204209    G4double sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT,
    … …  
    206211      //
    207212      // Cached trigonometric values
     213
     214    G4bool fPhiFullCone;
     215      //
     216      // Flag for identification of section or full cone
    208217};
    209218

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Cons.icc,v 1.8 2008/11/06 10:55:40 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Cons.icc,v 1.11 2009/06/09 16:08:23 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    3737
    3838inline
    39 void G4Cons::Initialise()
    40 {
    41   fCubicVolume= 0.;
    42   fSurfaceArea= 0.;
     39G4double G4Cons::GetInnerRadiusMinusZ() const
     40{
     41  return fRmin1 ;
     42}
     43
     44inline
     45G4double G4Cons::GetOuterRadiusMinusZ() const
     46{
     47  return fRmax1 ;
     48}
     49
     50inline
     51G4double G4Cons::GetInnerRadiusPlusZ() const
     52{
     53  return fRmin2 ;
     54}
     55
     56inline
     57G4double G4Cons::GetOuterRadiusPlusZ() const
     58{
     59  return fRmax2 ;
     60}
     61
     62inline
     63G4double G4Cons::GetZHalfLength() const
     64{
     65  return fDz ;
     66}
     67
     68inline 
     69G4double G4Cons::GetStartPhiAngle() const
     70{
     71  return fSPhi ;
     72}
     73
     74inline
     75G4double G4Cons::GetDeltaPhiAngle() const
     76{
     77  return fDPhi;
     78}
     79
     80inline
     81void G4Cons::Initialize()
     82{
     83  fCubicVolume = 0.;
     84  fSurfaceArea = 0.;
    4385  fpPolyhedron = 0;
    4486}
    … …  
    61103}
    62104
    63 inline
    64 G4double G4Cons::GetInnerRadiusMinusZ() const
    65 {
    66   return fRmin1 ;
    67 }
    68 
    69 inline
    70 G4double G4Cons::GetOuterRadiusMinusZ() const
    71 {
    72   return fRmax1 ;
    73 }
    74 
    75 inline
    76 G4double G4Cons::GetInnerRadiusPlusZ() const
    77 {
    78   return fRmin2 ;
    79 }
    80 
    81 inline
    82 G4double G4Cons::GetOuterRadiusPlusZ() const
    83 {
    84   return fRmax2 ;
    85 }
    86 
    87 inline
    88 G4double G4Cons::GetZHalfLength() const
    89 {
    90   return fDz ;
    91 }
    92 
    93 inline 
    94 G4double G4Cons::GetStartPhiAngle() const
    95 {
    96   return fSPhi ;
    97 }
    98 
    99 inline
    100 G4double G4Cons::GetDeltaPhiAngle() const
    101 {
    102   return fDPhi;
     105inline void G4Cons::CheckSPhiAngle(G4double sPhi)
     106{
     107  // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
     108
     109  if ( sPhi < 0 )
     110  {
     111    fSPhi = twopi - std::fmod(std::fabs(sPhi),twopi);
     112  }
     113  else
     114  {
     115    fSPhi = std::fmod(sPhi,twopi) ;
     116  }
     117  if ( fSPhi+fDPhi > twopi )
     118  {
     119    fSPhi -= twopi ;
     120  }
     121}
     122
     123inline void G4Cons::CheckDPhiAngle(G4double dPhi)
     124{
     125  fPhiFullCone = true;
     126  if ( dPhi >= twopi-kAngTolerance*0.5 )
     127  {
     128    fDPhi=twopi;
     129    fSPhi=0;
     130  }
     131  else
     132  {
     133    fPhiFullCone = false;
     134    if ( dPhi > 0 )
     135    {
     136      fDPhi = dPhi;
     137    }
     138    else
     139    {
     140      G4cerr << "ERROR - G4Cons()::CheckDPhiAngle()" << G4endl
     141             << "        Negative or zero delta-Phi (" << dPhi << ") in solid: "
     142             << GetName() << G4endl;
     143      G4Exception("G4Cons::CheckDPhiAngle()", "InvalidSetup",
     144                  FatalException, "Invalid dphi.");
     145    }
     146  }
     147}
     148
     149inline void G4Cons::CheckPhiAngles(G4double sPhi, G4double dPhi)
     150{
     151  CheckDPhiAngle(dPhi);
     152  if ( (fDPhi<twopi) && (sPhi) ) { CheckSPhiAngle(sPhi); }
     153  InitializeTrigonometry();
    103154}
    104155
    … …  
    107158{
    108159  fRmin1= Rmin1 ;
    109   Initialise();
     160  Initialize();
    110161}
    111162
    … …  
    114165{
    115166  fRmax1= Rmax1 ;
    116   Initialise();
     167  Initialize();
    117168}
    118169
    … …  
    121172{
    122173  fRmin2= Rmin2 ;
    123   Initialise();
     174  Initialize();
    124175}
    125176
    … …  
    128179{
    129180  fRmax2= Rmax2 ;
    130   Initialise();
     181  Initialize();
    131182}
    132183
    … …  
    135186{
    136187  fDz= newDz ;
    137   Initialise();
    138 }
    139 
    140 inline
    141 void G4Cons::SetStartPhiAngle ( G4double newSPhi )
    142 {
    143   fSPhi= newSPhi;
    144   Initialise();
    145   InitializeTrigonometry();
     188  Initialize();
     189}
     190
     191inline
     192void G4Cons::SetStartPhiAngle ( G4double newSPhi, G4bool compute )
     193{
     194  // Flag 'compute' can be used to explicitely avoid recomputation of
     195  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
     196
     197  CheckSPhiAngle(newSPhi);
     198  fPhiFullCone = false;
     199  if (compute)  { InitializeTrigonometry(); }
     200  Initialize();
    146201}
    147202
    148203void G4Cons::SetDeltaPhiAngle ( G4double newDPhi )
    149204{
    150   if ( newDPhi >= twopi-kAngTolerance*0.5 )
    151   {
    152     fPhiFullCone = true;
    153   }
    154   else if ( newDPhi > 0 )
    155   {
    156     fPhiFullCone = false;
    157   }
    158   else
    159   {
    160      G4cerr << "ERROR - G4Cons()::SetDeltaPhiAngle() : " << GetName() << G4endl
    161             << "        Negative delta-Phi ! - " << newDPhi << G4endl;
    162      G4Exception("G4Cons::SetDeltaPhiAngle()", "InvalidSetup",
    163                   FatalException, "Invalid dphi.");
    164   }
    165   fDPhi=  newDPhi;
    166   Initialise();
    167   InitializeTrigonometry();
     205  CheckPhiAngles(fSPhi, newDPhi);
     206  Initialize();
    168207}
    169208

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

    r1058 r1228  
    2626//
    2727// $Id: G4Orb.hh,v 1.11 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4Orb.icc,v 1.5 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2626//
    2727// $Id: G4Para.hh,v 1.18 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4Para.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Sphere.hh,v 1.21 2008/11/21 09:50:05 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Sphere.hh,v 1.24 2009/03/31 07:51:49 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    3636// Class description:
    3737//
    38 //   A G4Sphere is, in the general case, section of a spherical shell,
     38//   A G4Sphere is, in the general case, a section of a spherical shell,
    3939//   between specified phi and theta angles
    4040//
    … …  
    8383                   G4double pSPhi, G4double pDPhi,
    8484                   G4double pSTheta, G4double pDTheta);
     85      //
     86      // Constructs a sphere or sphere shell section
     87      // with the given name and dimensions
    8588       
    86     virtual ~G4Sphere() ;
    87    
     89   ~G4Sphere();
     90      //
     91      // Destructor
     92
    8893    // Accessors
    8994       
    … …  
    99104    inline void SetInnerRadius    (G4double newRMin);
    100105    inline void SetOuterRadius    (G4double newRmax);
    101     inline void SetStartPhiAngle  (G4double newSphi);
     106    inline void SetStartPhiAngle  (G4double newSphi, G4bool trig=true);
    102107    inline void SetDeltaPhiAngle  (G4double newDphi);
    103108    inline void SetStartThetaAngle(G4double newSTheta);
    … …  
    107112
    108113    inline G4double GetCubicVolume();
    109     inline G4double GetSurfaceArea();
     114    G4double GetSurfaceArea();
    110115
    111116    void ComputeDimensions(      G4VPVParameterisation* p,
    … …  
    142147
    143148    // Visualisation functions
     149
    144150    G4VisExtent   GetExtent          () const;   
    145151    void          DescribeYourselfTo(G4VGraphicsScene& scene) const;
    … …  
    150156   
    151157    G4Sphere(__void__&);
     158      //
    152159      // Fake default constructor for usage restricted to direct object
    153160      // persistency for clients requiring preallocation of memory for
    … …  
    165172    inline void SetInsideRadius(G4double newRmin);
    166173
    167   protected:
     174  private:
    168175 
    169176    G4ThreeVectorList*
    170177    CreateRotatedVertices(const G4AffineTransform& pTransform,
    171178                                G4int& noPolygonVertices) const;
    172  
     179      //
     180      // Creates the List of transformed vertices in the format required
     181      // for G4VSolid:: ClipCrossSection and ClipBetweenSections
     182
     183    inline void Initialize();
     184      //
     185      // Reset relevant values to zero
     186
     187    inline void CheckThetaAngles(G4double sTheta, G4double dTheta);
     188    inline void CheckSPhiAngle(G4double sPhi);
     189    inline void CheckDPhiAngle(G4double dPhi);
     190    inline void CheckPhiAngles(G4double sPhi, G4double dPhi);
     191      //
     192      // Reset relevant flags and angle values
     193
     194    inline void InitializePhiTrigonometry();
     195    inline void InitializeThetaTrigonometry();
     196      //
     197      // Recompute relevant trigonometric values and cache them
     198
     199    G4ThreeVector ApproxSurfaceNormal( const G4ThreeVector& p) const;
     200      //
     201      // Algorithm for SurfaceNormal() following the original
     202      // specification for points not on the surface
     203
     204  private:
     205
    173206    // Used by distanceToOut
    174  
     207    //
    175208    enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta};
    176209 
    177210    // used by normal
    178  
     211    //
    179212    enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};
    180213
    181   private:
    182 
    183     G4ThreeVector ApproxSurfaceNormal( const G4ThreeVector& p) const;
    184       // Algorithm for SurfaceNormal() following the original
    185       // specification for points not on the surface
    186 
    187   private:
    188 
    189     G4double kAngTolerance, kRadTolerance;
    190 
    191     G4double fRmin,fRmax,
    192              fSPhi,fDPhi,
    193              fSTheta,fDTheta;
    194     G4double fEpsilon;
     214    G4double fEpsilon, fRminTolerance, fRmaxTolerance, kAngTolerance;
     215      //
     216      // Radial and angular tolerances
     217
     218    G4double fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta;
     219      //
     220      // Radial and angular dimensions
     221
     222    G4double sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT,
     223             sinSPhi, cosSPhi, sinEPhi, cosEPhi, hDPhi, cPhi, ePhi;
     224      //
     225      // Cached trigonometric values for Phi angle
     226
     227    G4double sinSTheta, cosSTheta, sinETheta, cosETheta,
     228             tanSTheta, tanSTheta2, tanETheta, tanETheta2, eTheta;
     229      //
     230      // Cached trigonometric values for Theta angle
     231
     232    G4bool fFullPhiSphere, fFullThetaSphere, fFullSphere;
     233      //
     234      // Flags for identification of section, shell or full sphere
    195235};
    196236

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Sphere.icc,v 1.8 2008/11/21 09:50:05 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Sphere.icc,v 1.11 2009/05/13 11:23:00 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    7777}
    7878
     79inline
     80void G4Sphere::Initialize()
     81{
     82  fCubicVolume = 0.;
     83  fSurfaceArea = 0.;
     84  fpPolyhedron = 0;
     85}
     86
     87inline
     88void G4Sphere::InitializePhiTrigonometry()
     89{
     90  hDPhi = 0.5*fDPhi;                       // half delta phi
     91  cPhi  = fSPhi + hDPhi;
     92  ePhi  = fSPhi + fDPhi;
     93
     94  sinCPhi    = std::sin(cPhi);
     95  cosCPhi    = std::cos(cPhi);
     96  cosHDPhiIT = std::cos(hDPhi - 0.5*kAngTolerance); // inner/outer tol half dphi
     97  cosHDPhiOT = std::cos(hDPhi + 0.5*kAngTolerance);
     98  sinSPhi = std::sin(fSPhi);
     99  cosSPhi = std::cos(fSPhi);
     100  sinEPhi = std::sin(ePhi);
     101  cosEPhi = std::cos(ePhi);
     102}
     103
     104inline
     105void G4Sphere::InitializeThetaTrigonometry()
     106{
     107  eTheta  = fSTheta + fDTheta;
     108
     109  sinSTheta = std::sin(fSTheta);
     110  cosSTheta = std::cos(fSTheta);
     111  sinETheta = std::sin(eTheta);
     112  cosETheta = std::cos(eTheta);
     113
     114  tanSTheta  = std::tan(fSTheta);
     115  tanSTheta2 = tanSTheta*tanSTheta;
     116  tanETheta  = std::tan(eTheta);
     117  tanETheta2 = tanETheta*tanETheta;
     118}
     119
     120inline
     121void G4Sphere::CheckThetaAngles(G4double sTheta, G4double dTheta)
     122{
     123  if ( (sTheta<0) || (sTheta>pi) )
     124  {
     125    G4cerr << "ERROR - G4Sphere()::CheckThetaAngles()" << G4endl
     126           << "        Invalid starting Theta angle for solid: " << GetName()
     127           << G4endl;
     128    G4Exception("G4Sphere::CheckThetaAngles()", "InvalidSetup",
     129                FatalException, "sTheta outside 0-PI range.");
     130  }
     131  else
     132  {
     133    fSTheta=sTheta;
     134  }
     135  if ( dTheta+sTheta >= pi )
     136  {
     137    fDTheta=pi-sTheta;
     138  }
     139  else if ( dTheta > 0 )
     140  {
     141    fDTheta=dTheta;
     142  }
     143  else
     144  {
     145    G4cerr << "ERROR - G4Sphere()::CheckThetaAngles(): " << G4endl
     146           << "        Negative delta-Theta (" << dTheta << "), for solid: "
     147           << GetName() << G4endl;
     148    G4Exception("G4Sphere::CheckThetaAngles()", "InvalidSetup",
     149                FatalException, "Invalid dTheta.");
     150  }
     151  if ( fDTheta-fSTheta < pi ) { fFullThetaSphere = false; }
     152  else                        { fFullThetaSphere = true ; }
     153  fFullSphere = fFullPhiSphere && fFullThetaSphere;
     154
     155  InitializeThetaTrigonometry();
     156}
     157
     158inline
     159void G4Sphere::CheckSPhiAngle(G4double sPhi)
     160{
     161  // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
     162
     163  if ( sPhi < 0 )
     164  {
     165    fSPhi = twopi - std::fmod(std::fabs(sPhi),twopi);
     166  }
     167  else
     168  {
     169    fSPhi = std::fmod(sPhi,twopi) ;
     170  }
     171  if ( fSPhi+fDPhi > twopi )
     172  {
     173    fSPhi -= twopi ;
     174  }
     175}
     176
     177inline
     178void G4Sphere::CheckDPhiAngle(G4double dPhi)
     179{
     180  fFullPhiSphere = true;
     181  if ( dPhi >= twopi-kAngTolerance*0.5 )
     182  {
     183    fDPhi=twopi;
     184    fSPhi=0;
     185  }
     186  else
     187  {
     188    fFullPhiSphere = false;
     189    if ( dPhi > 0 )
     190    {
     191      fDPhi = dPhi;
     192    }
     193    else
     194    {
     195      G4cerr << "ERROR - G4Sphere()::CheckDPhiAngle(): " << GetName() << G4endl
     196             << "        Negative delta-Phi ! - "
     197             << dPhi << G4endl;
     198      G4Exception("G4Sphere::CheckDPhiAngle()", "InvalidSetup",
     199                  FatalException, "Invalid dphi.");
     200    }
     201  }
     202}
     203
     204inline
     205void G4Sphere::CheckPhiAngles(G4double sPhi, G4double dPhi)
     206{
     207  CheckDPhiAngle(dPhi);
     208  if (!fFullPhiSphere && sPhi) { CheckSPhiAngle(sPhi); }
     209  fFullSphere = fFullPhiSphere && fFullThetaSphere;
     210
     211  InitializePhiTrigonometry();
     212}
     213
    79214inline
    80215void G4Sphere::SetInsideRadius(G4double newRmin)
    81216{
    82217  fRmin= newRmin;
    83   fCubicVolume= 0.;
    84   fSurfaceArea= 0.;
    85   fpPolyhedron = 0;
     218  Initialize();
    86219}
    87220
    … …  
    90223{
    91224  fRmin= newRmin;
    92   fCubicVolume= 0.;
    93   fSurfaceArea= 0.;
    94   fpPolyhedron = 0;
     225  Initialize();
    95226}
    96227
    … …  
    99230{
    100231  fRmax= newRmax;
    101   fCubicVolume= 0.;
    102   fSurfaceArea= 0.;
    103   fpPolyhedron = 0;
    104 }
    105 
    106 inline
    107 void G4Sphere::SetStartPhiAngle(G4double newSphi)
    108 {
    109   fSPhi= newSphi;
    110   fCubicVolume= 0.;
    111   fSurfaceArea= 0.;
    112   fpPolyhedron = 0;
    113 }
    114 
    115 inline
    116 void G4Sphere::SetDeltaPhiAngle(G4double newDphi)
    117 {
    118   fDPhi= newDphi;
    119   fCubicVolume= 0.;
    120   fSurfaceArea= 0.;
    121   fpPolyhedron = 0;
     232  Initialize();
     233}
     234
     235inline
     236void G4Sphere::SetStartPhiAngle(G4double newSPhi, G4bool compute)
     237{
     238  // Flag 'compute' can be used to explicitely avoid recomputation of
     239  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
     240
     241  CheckSPhiAngle(newSPhi);
     242  fFullPhiSphere = false;
     243  if (compute)  { InitializePhiTrigonometry(); }
     244  Initialize();
     245}
     246
     247inline
     248void G4Sphere::SetDeltaPhiAngle(G4double newDPhi)
     249{
     250  CheckPhiAngles(fSPhi, newDPhi);
     251  Initialize();
    122252}
    123253
    … …  
    125255void G4Sphere::SetStartThetaAngle(G4double newSTheta)
    126256{
    127   fSTheta=newSTheta;
    128   fCubicVolume= 0.;
    129   fSurfaceArea= 0.;
    130   fpPolyhedron = 0;
     257  CheckThetaAngles(newSTheta, fDTheta);
     258  Initialize();
    131259}
    132260
    … …  
    134262void G4Sphere::SetDeltaThetaAngle(G4double newDTheta)
    135263{
    136   fDTheta=newDTheta;
    137   fCubicVolume= 0.;
    138   fSurfaceArea= 0.;
    139   fpPolyhedron = 0;
     264  CheckThetaAngles(fSTheta, newDTheta);
     265  Initialize();
    140266}
    141267
    … …  
    182308{
    183309  if(fCubicVolume != 0.) {;}
    184   else   fCubicVolume = fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta))*
    185                               (fRmax*fRmax*fRmax-fRmin*fRmin*fRmin)/3.;
     310  else { fCubicVolume = fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta))*
     311                              (fRmax*fRmax*fRmax-fRmin*fRmin*fRmin)/3.; }
    186312  return fCubicVolume;
    187313}
    188 
    189 
    190 inline
    191 G4double G4Sphere::GetSurfaceArea()
    192 {
    193   if(fSurfaceArea != 0.) {;}
    194   else
    195   {   
    196     G4double Rsq=fRmax*fRmax;
    197     G4double rsq=fRmin*fRmin;
    198          
    199     fSurfaceArea = fDPhi*(rsq+Rsq)*(std::cos(fSTheta)
    200                  - std::cos(fSTheta+fDTheta));
    201     if(fDPhi < twopi)
    202     {
    203       fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
    204     }
    205     if(fSTheta >0)
    206     {
    207       G4double acos1=std::acos( std::pow(std::sin(fSTheta),2)
    208                                *std::cos(fDPhi)
    209                                +std::pow(std::cos(fSTheta),2));
    210       if(fDPhi>pi)
    211       {
    212         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
    213       }
    214       else
    215       {
    216         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
    217       }
    218     }
    219     if(fDTheta+fSTheta < pi)
    220     {
    221       G4double acos2=std::acos( std::pow(std::sin(fSTheta+fDTheta),2)
    222                                *std::cos(fDPhi)
    223                                +std::pow(std::cos(fSTheta+fDTheta),2));
    224       if(fDPhi>pi)
    225       {
    226         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
    227       }
    228       else
    229       {
    230         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
    231       }
    232     }
    233   }
    234   return fSurfaceArea;
    235 }

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

    r1058 r1228  
    2626//
    2727// $Id: G4Torus.hh,v 1.27 2007/05/18 07:38:00 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4Torus.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trap.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trap.icc,v 1.8 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trd.hh,v 1.16 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trd.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Tubs.hh,v 1.21 2008/11/06 10:55:40 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Tubs.hh,v 1.22 2009/03/26 16:25:44 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    8686      // Constructs a tubs with the given name and dimensions
    8787
    88     virtual ~G4Tubs();
     88   ~G4Tubs();
    8989      //
    9090      // Destructor
    … …  
    103103    inline void SetOuterRadius   (G4double newRMax);
    104104    inline void SetZHalfLength   (G4double newDz);
    105     inline void SetStartPhiAngle (G4double newSPhi);
     105    inline void SetStartPhiAngle (G4double newSPhi, G4bool trig=true);
    106106    inline void SetDeltaPhiAngle (G4double newDPhi);
    107107   
    … …  
    159159    inline G4double GetDPhi() const;
    160160
    161   protected:
     161  private:
    162162
    163163    G4ThreeVectorList*
    … …  
    167167      // for G4VSolid:: ClipCrossSection and ClipBetweenSections
    168168
    169     G4double fRMin, fRMax, fDz, fSPhi, fDPhi;
    170     G4bool fPhiFullTube;
    171    
    172       // Used by distanceToOut
    173 
    174     enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
    175 
    176       // Used by normal
    177 
    178     enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
    179 
    180   private:
    181 
    182169    inline void Initialize();
    183170      //
    184171      // Reset relevant values to zero
     172
     173    inline void CheckSPhiAngle(G4double sPhi);
     174    inline void CheckDPhiAngle(G4double dPhi);
     175    inline void CheckPhiAngles(G4double sPhi, G4double dPhi);
     176      //
     177      // Reset relevant flags and angle values
    185178
    186179    inline void InitializeTrigonometry();
    … …  
    195188  private:
    196189
     190    // Used by distanceToOut
     191    //
     192    enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
     193
     194    // Used by normal
     195    //
     196    enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
     197
    197198    G4double kRadTolerance, kAngTolerance;
    198199      //
    199200      // Radial and angular tolerances
    200201
     202    G4double fRMin, fRMax, fDz, fSPhi, fDPhi;
     203      //
     204      // Radial and angular dimensions
     205   
    201206    G4double sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT,
    202207             sinSPhi, cosSPhi, sinEPhi, cosEPhi;
    203208      //
    204209      // Cached trigonometric values
     210
     211    G4bool fPhiFullTube;
     212      //
     213      // Flag for identification of section or full tube
    205214};
    206215

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Tubs.icc,v 1.11 2008/11/06 10:55:40 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Tubs.icc,v 1.14 2009/06/09 16:08:23 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    9191}
    9292
     93inline void G4Tubs::CheckSPhiAngle(G4double sPhi)
     94{
     95  // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
     96
     97  if ( sPhi < 0 )
     98  {
     99    fSPhi = twopi - std::fmod(std::fabs(sPhi),twopi);
     100  }
     101  else
     102  {
     103    fSPhi = std::fmod(sPhi,twopi) ;
     104  }
     105  if ( fSPhi+fDPhi > twopi )
     106  {
     107    fSPhi -= twopi ;
     108  }
     109}
     110
     111inline void G4Tubs::CheckDPhiAngle(G4double dPhi)
     112{
     113  fPhiFullTube = true;
     114  if ( dPhi >= twopi-kAngTolerance*0.5 )
     115  {
     116    fDPhi=twopi;
     117    fSPhi=0;
     118  }
     119  else
     120  {
     121    fPhiFullTube = false;
     122    if ( dPhi > 0 )
     123    {
     124      fDPhi = dPhi;
     125    }
     126    else
     127    {
     128      G4cerr << "ERROR - G4Tubs()::CheckDPhiAngle()" << G4endl
     129             << "        Negative or zero delta-Phi (" << dPhi << ") in solid: "
     130             << GetName() << G4endl;
     131      G4Exception("G4Tubs::CheckDPhiAngle()", "InvalidSetup",
     132                  FatalException, "Invalid dphi.");
     133    }
     134  }
     135}
     136
     137inline void G4Tubs::CheckPhiAngles(G4double sPhi, G4double dPhi)
     138{
     139  CheckDPhiAngle(dPhi);
     140  if ( (fDPhi<twopi) && (sPhi) ) { CheckSPhiAngle(sPhi); }
     141  InitializeTrigonometry();
     142}
     143
    93144inline
    94145void G4Tubs::SetInnerRadius (G4double newRMin)
    … …  
    113164
    114165inline
    115 void G4Tubs::SetStartPhiAngle (G4double newSPhi)
    116 {
    117   fSPhi= newSPhi;
    118   Initialize();
    119   InitializeTrigonometry();
     166void G4Tubs::SetStartPhiAngle (G4double newSPhi, G4bool compute)
     167{
     168  // Flag 'compute' can be used to explicitely avoid recomputation of
     169  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
     170
     171  CheckSPhiAngle(newSPhi);
     172  fPhiFullTube = false;
     173  if (compute)  { InitializeTrigonometry(); }
     174  Initialize();
    120175}
    121176
    … …  
    123178void G4Tubs::SetDeltaPhiAngle (G4double newDPhi)
    124179{
    125   if ( newDPhi >= twopi-kAngTolerance*0.5 )
    126   {
    127     fPhiFullTube = true;
    128   }
    129   else if ( newDPhi > 0 )
    130   {
    131     fPhiFullTube = false;
    132   }
    133   else
    134   {
    135      G4cerr << "ERROR - G4Tubs()::SetDeltaPhiAngle() : " << GetName() << G4endl
    136             << "        Negative delta-Phi ! - " << newDPhi << G4endl;
    137      G4Exception("G4Tubs::SetDeltaPhiAngle()", "InvalidSetup",
    138                   FatalException, "Invalid dphi.");
    139   }
    140   fDPhi= newDPhi;
    141   Initialize();
    142   InitializeTrigonometry();
     180  CheckPhiAngles(fSPhi, newDPhi);
     181  Initialize();
    143182}
    144183

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

    r1058 r1228  
    2626//
    2727// $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2626//
    2727// $Id: G4CSGSolid.cc,v 1.13 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// --------------------------------------------------------------------

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Cons.cc,v 1.60 2008/11/06 15:26:53 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Cons.cc,v 1.67 2009/11/12 11:53:11 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    3535// History:
    3636//
     37// 12.10.09 T.Nikitina: Added to DistanceToIn(p,v) check on the direction in
     38//                      case of point on surface
    3739// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
    3840// 13.09.96 V.Grichine: Review and final modifications
    … …  
    7981                      G4double pDz,
    8082                      G4double pSPhi, G4double pDPhi)
    81   : G4CSGSolid(pName)
     83  : G4CSGSolid(pName), fSPhi(0), fDPhi(0)
    8284{
    83   // Check z-len
    84 
    8585  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
    8686  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
    8787
     88  // Check z-len
     89  //
    8890  if ( pDz > 0 )
    8991  {
    … …  
    100102
    101103  // Check radii
    102 
     104  //
    103105  if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) )
    104106  {
    … …  
    121123  }
    122124
    123   fPhiFullCone = true;
    124   if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
    125   {
    126     fDPhi=twopi;
    127     fSPhi=0;
    128   }
    129   else
    130   {
    131     fPhiFullCone = false;
    132     if ( pDPhi > 0 )
    133     {
    134       fDPhi = pDPhi;
    135     }
    136     else
    137     {
    138       G4cerr << "ERROR - G4Cons()::G4Cons(): " << GetName() << G4endl
    139              << "        Negative delta-Phi ! - "
    140              << pDPhi << G4endl;
    141       G4Exception("G4Cons::G4Cons()", "InvalidSetup",
    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();
     125  // Check angles
     126  //
     127  CheckPhiAngles(pSPhi, pDPhi);
    161128}
    162129
    … …  
    705672{
    706673  G4double snxt = kInfinity ;      // snxt = default return value
    707 
     674  const G4double dRmax = 100*std::min(fRmax1,fRmax2);
    708675  static const G4double halfCarTolerance=kCarTolerance*0.5;
    709676  static const G4double halfRadTolerance=kRadTolerance*0.5;
    … …  
    717684  G4double tolODz,tolIDz ;
    718685
    719   G4double Dist,s,xi,yi,zi,ri=0.,rhoi2,cosPsi ; // Intersection point variables
     686  G4double Dist,s,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
    720687
    721688  G4double t1,t2,t3,b,c,d ;    // Quadratic solver variables
    722689  G4double nt1,nt2,nt3 ;
    723690  G4double Comp ;
     691
     692  G4ThreeVector Normal;
    724693
    725694  // Cone Precalcs
    … …  
    887856        if ( s > 0 )  // If 'forwards'. Check z intersection
    888857        {
     858          if ( s>dRmax ) // Avoid rounding errors due to precision issues on
     859          {              // 64 bits systems. Split long distances and recompute
     860            G4double fTerm = s-std::fmod(s,dRmax);
     861            s = fTerm + DistanceToIn(p+fTerm*v,v);
     862          }
    889863          zi = p.z() + s*v.z() ;
    890864
    … …  
    917891      {
    918892        // Inside cones, delta r -ve, inside z extent
    919 
     893        // Point is on the Surface => check Direction using  Normal.dot(v)
     894
     895        xi     = p.x() ;
     896        yi     = p.y()  ;
     897        risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
     898        Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
    920899        if ( !fPhiFullCone )
    921900        {
    922901          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
    923 
    924           if (cosPsi >= cosHDPhiIT)  { return 0.0; }
    925         }
    926         else  { return 0.0; }
     902          if ( cosPsi >= cosHDPhiIT )
     903          {
     904            if ( Normal.dot(v) <= 0 )  { return 0.0; }
     905          }
     906        }
     907        else
     908        {             
     909          if ( Normal.dot(v) <= 0 )  { return 0.0; }
     910        }
    927911      }
    928912    }
    … …  
    972956
    973957  if (rMinAv)
    974   {
     958  { 
    975959    nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
    976960    nt2 = t2 - tanRMin*v.z()*rin ;
    … …  
    993977          if ( s >= 0 )   // > 0
    994978          {
     979            if ( s>dRmax ) // Avoid rounding errors due to precision issues on
     980            {              // 64 bits systems. Split long distance and recompute
     981              G4double fTerm = s-std::fmod(s,dRmax);
     982              s = fTerm + DistanceToIn(p+fTerm*v,v);
     983            }
    995984            zi = p.z() + s*v.z() ;
    996985
    … …  
    1004993                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    1005994
    1006                 if (cosPsi >= cosHDPhiIT)  { snxt = s; }
     995                if (cosPsi >= cosHDPhiIT)
     996                {
     997                  if ( s > halfRadTolerance )  { snxt=s; }
     998                  else
     999                  {
     1000                    // Calculate a normal vector in order to check Direction
     1001
     1002                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1003                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
     1004                    if ( Normal.dot(v) <= 0 )  { snxt = s; }
     1005                  }
     1006                }
    10071007              }
    1008               else  { return s; }
     1008              else
     1009              {
     1010                if ( s > halfRadTolerance )  { return s; }
     1011                else
     1012                {
     1013                  // Calculate a normal vector in order to check Direction
     1014
     1015                  xi     = p.x() + s*v.x() ;
     1016                  yi     = p.y() + s*v.y() ;
     1017                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1018                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
     1019                  if ( Normal.dot(v) <= 0 )  { return s; }
     1020                }
     1021              }
    10091022            }
    10101023          }
    … …  
    10321045            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s > 0
    10331046            {
     1047              if ( s>dRmax ) // Avoid rounding errors due to precision issues
     1048              {              // seen on 64 bits systems. Split and recompute
     1049                G4double fTerm = s-std::fmod(s,dRmax);
     1050                s = fTerm + DistanceToIn(p+fTerm*v,v);
     1051              }
    10341052              if ( !fPhiFullCone )
    10351053              {
    … …  
    10381056                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    10391057
    1040                 if (cosPsi >= cosHDPhiOT)  { snxt = s; }
     1058                if (cosPsi >= cosHDPhiOT)
     1059                {
     1060                  if ( s > halfRadTolerance )  { snxt=s; }
     1061                  else
     1062                  {
     1063                    // Calculate a normal vector in order to check Direction
     1064
     1065                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1066                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
     1067                    if ( Normal.dot(v) <= 0 )  { snxt = s; }
     1068                  }
     1069                }
    10411070              }
    1042               else  { return s; }
     1071              else
     1072              {
     1073                if( s > halfRadTolerance )  { return s; }
     1074                else
     1075                {
     1076                  // Calculate a normal vector in order to check Direction
     1077
     1078                  xi     = p.x() + s*v.x() ;
     1079                  yi     = p.y() + s*v.y() ;
     1080                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1081                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
     1082                  if ( Normal.dot(v) <= 0 )  { return s; }
     1083                }
     1084              }
    10431085            }
    10441086          }
    … …  
    10511093            if ( (s >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // s>0
    10521094            {
     1095              if ( s>dRmax ) // Avoid rounding errors due to precision issues
     1096              {              // seen on 64 bits systems. Split and recompute
     1097                G4double fTerm = s-std::fmod(s,dRmax);
     1098                s = fTerm + DistanceToIn(p+fTerm*v,v);
     1099              }
    10531100              if ( !fPhiFullCone )
    10541101              {
    … …  
    10571104                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
    10581105
    1059                 if (cosPsi >= cosHDPhiIT)  { snxt = s; }
     1106                if (cosPsi >= cosHDPhiIT)
     1107                {
     1108                  if ( s > halfRadTolerance )  { snxt=s; }
     1109                  else
     1110                  {
     1111                    // Calculate a normal vector in order to check Direction
     1112
     1113                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1114                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
     1115                    if ( Normal.dot(v) <= 0 )  { snxt = s; }
     1116                  }
     1117                }
    10601118              }
    1061               else  { return s; }
     1119              else
     1120              {
     1121                if ( s > halfRadTolerance )  { return s; }
     1122                else
     1123                {
     1124                  // Calculate a normal vector in order to check Direction
     1125
     1126                  xi     = p.x() + s*v.x() ;
     1127                  yi     = p.y() + s*v.y() ;
     1128                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1129                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
     1130                  if ( Normal.dot(v) <= 0 )  { return s; }
     1131                }
     1132              }
    10621133            }
    10631134          }
    … …  
    10991170              zi = p.z() + s*v.z() ;
    11001171              ri = rMinAv + zi*tanRMin ;
    1101 
     1172             
    11021173              if ( ri > 0 )   // 2nd root
    11031174              {
    … …  
    11071178                if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
    11081179                {
     1180                  if ( s>dRmax ) // Avoid rounding errors due to precision issue
     1181                  {              // seen on 64 bits systems. Split and recompute
     1182                    G4double fTerm = s-std::fmod(s,dRmax);
     1183                    s = fTerm + DistanceToIn(p+fTerm*v,v);
     1184                  }
    11091185                  if ( !fPhiFullCone )
    11101186                  {
    … …  
    11361212            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
    11371213            {
     1214              if ( s>dRmax ) // Avoid rounding errors due to precision issues
     1215              {              // seen on 64 bits systems. Split and recompute
     1216                G4double fTerm = s-std::fmod(s,dRmax);
     1217                s = fTerm + DistanceToIn(p+fTerm*v,v);
     1218              }
    11381219              if ( !fPhiFullCone )
    11391220              {
    … …  
    13351416  // Vars for intersection within tolerance
    13361417
    1337   ESide    sidetol ;
     1418  ESide    sidetol = kNull ;
    13381419  G4double slentol = kInfinity ;
    13391420
    … …  
    16691750            xi     = p.x() + slentol*v.x() ;
    16701751            yi     = p.y() + slentol*v.y() ;
    1671             risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
    1672             Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMin/secRMin) ;
    1673            
     1752            if( sidetol==kRMax )
     1753            {
     1754              risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
     1755              Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
     1756            }
     1757            else
     1758            {
     1759              risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
     1760              Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
     1761            }
    16741762            if( Normal.dot(v) > 0 )
    16751763            {

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

    r1058 r1228  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Orb.cc,v 1.24 2007/05/18 07:38:01 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     26// $Id: G4Orb.cc,v 1.30 2009/11/30 10:20:38 gcosmo Exp $
     27// GEANT4 tag $Name: geant4-09-03 $
    2828//
    2929// class G4Orb
    … …  
    298298  rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z() ;
    299299
     300  G4double rad = std::sqrt(rad2);
     301
    300302  // G4double rad = std::sqrt(rad2);
    301303  // Check radial surface
    … …  
    304306  tolRMax = fRmax - fRmaxTolerance*0.5 ;
    305307   
    306   if ( rad2 <= tolRMax*tolRMax )  in = kInside ;
     308  if ( rad <= tolRMax )  { in = kInside ; }
    307309  else
    308310  {
    309311    tolRMax = fRmax + fRmaxTolerance*0.5 ;       
    310     if ( rad2 <= tolRMax*tolRMax ) in = kSurface ;
    311     else                           in = kOutside ;
     312    if ( rad <= tolRMax )  { in = kSurface ; }
     313    else                   { in = kOutside ; }
    312314  }
    313315  return in;
    … …  
    360362  G4double c, d2, s = kInfinity ;
    361363
     364  const G4double dRmax = 100.*fRmax;
     365
    362366  // General Precalcs
    363367
    … …  
    384388  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
    385389
    386   c = rad2 - fRmax*fRmax ;
     390
     391  G4double rad = std::sqrt(rad2);
     392  c = (rad - fRmax)*(rad + fRmax);
    387393
    388394  if ( c > fRmaxTolerance*fRmax )
    … …  
    396402    {
    397403      s = -pDotV3d - std::sqrt(d2) ;
    398 
    399       if (s >= 0 ) return snxt = s;
    400            
     404      if ( s >= 0 )
     405      {
     406        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
     407        {              // 64 bits systems. Split long distances and recompute
     408          G4double fTerm = s-std::fmod(s,dRmax);
     409          s = fTerm + DistanceToIn(p+fTerm*v,v);
     410        }
     411        return snxt = s;
     412      }
    401413    }
    402414    else    // No intersection with G4Orb
    … …  
    410422    {
    411423      d2 = pDotV3d*pDotV3d - c ;             
    412       //  if ( pDotV3d >= 0 ) return snxt = kInfinity;
    413       if ( d2 < fRmaxTolerance*fRmax || pDotV3d >= 0 ) return snxt = kInfinity;
    414       else                return snxt = 0.;
    415     }
     424      if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
     425      {
     426        return snxt = kInfinity;
     427      }
     428      else
     429      {
     430        return snxt = 0.;
     431      }
     432    }
     433#ifdef G4CSGDEBUG
    416434    else // inside ???
    417435    {
    … …  
    419437                  JustWarning, "Point p is inside !?");
    420438    }
     439#endif
    421440  }
    422441  return snxt;
    … …  
    431450G4double G4Orb::DistanceToIn( const G4ThreeVector& p ) const
    432451{
    433   G4double safe=0.0, rad  = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
    434                  safe = rad - fRmax;
    435   if( safe < 0 ) safe = 0. ;
     452  G4double safe = 0.0,
     453           rad  = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
     454  safe = rad - fRmax;
     455  if( safe < 0 ) { safe = 0.; }
    436456  return safe;
    437457}
    … …  
    475495 
    476496  const G4double  Rmax_plus = fRmax + fRmaxTolerance*0.5;
    477 
    478   if( rad2 <= Rmax_plus*Rmax_plus )
    479   {
    480     c = rad2-fRmax*fRmax ;
    481 
    482     if ( c < fRmaxTolerance*fRmax)
     497  G4double rad = std::sqrt(rad2);
     498
     499  if ( rad <= Rmax_plus )
     500  {
     501    c = (rad - fRmax)*(rad + fRmax);
     502
     503    if ( c < fRmaxTolerance*fRmax )
    483504    {
    484505      // Within tolerant Outer radius

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

    r1058 r1228  
    2626//
    2727// $Id: G4Para.cc,v 1.39 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// class G4Para

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Sphere.cc,v 1.68 2008/07/07 09:35:16 grichine Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Sphere.cc,v 1.84 2009/08/07 15:56:23 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// class G4Sphere
    … …  
    3434// History:
    3535//
     36// 14.09.09 T.Nikitina: fix for phi section in DistanceToOut(p,v,..),as for G4Tubs,G4Cons
     37// 26.03.09 G.Cosmo   : optimisations and uniform use of local radial tolerance
    3638// 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...)
    3739// 22.07.05 O.Link    : Added check for intersection with double cone
    … …  
    4143// 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside
    4244// 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections
    43 // 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi < SPhi, SPhi<0
     45// 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi<SPhi, SPhi<0
    4446// 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
    4547// 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
    … …  
    5355// --------------------------------------------------------------------
    5456
    55 #include <assert.h>
    56 
    5757#include "G4Sphere.hh"
    5858
    … …  
    9292                          G4double pSPhi, G4double pDPhi,
    9393                          G4double pSTheta, G4double pDTheta )
    94   : G4CSGSolid(pName)
     94  : G4CSGSolid(pName), fFullPhiSphere(true), fFullThetaSphere(true)
    9595{
    96   fEpsilon = 1.0e-14;
    97 
    98   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
     96  fEpsilon = 2.0e-11;  // relative radial tolerance constant
     97
    9998  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
    10099
    101   // Check radii
    102 
    103   if (pRmin<pRmax&&pRmin>=0)
     100  // Check radii and set radial tolerances
     101
     102  G4double kRadTolerance = G4GeometryTolerance::GetInstance()
     103                         ->GetRadialTolerance();
     104  if ( (pRmin < pRmax) && (pRmax >= 10*kRadTolerance) && (pRmin >= 0) )
    104105  {
    105106    fRmin=pRmin; fRmax=pRmax;
     107    fRminTolerance = (pRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
     108    fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
    106109  }
    107110  else
    … …  
    116119  // Check angles
    117120
    118   if (pDPhi>=twopi)
    119   {
    120     fDPhi=twopi;
    121   }
    122   else if (pDPhi>0)
    123   {
    124     fDPhi=pDPhi;
    125   }
    126   else
    127   {
    128     G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
    129            << "        Negative Z delta-Phi ! - "
    130            << pDPhi << G4endl;
    131     G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
    132                 "Invalid DPhi.");
    133   }
    134 
    135   // Convert fSPhi to 0-2PI
    136 
    137   if (pSPhi<0)
    138   {
    139     fSPhi=twopi-std::fmod(std::fabs(pSPhi),twopi);
    140   }
    141   else
    142   {
    143     fSPhi=std::fmod(pSPhi,twopi);
    144   }
    145 
    146   // Sphere is placed such that fSPhi+fDPhi>twopi !
    147   // fSPhi could be < 0 !!?
    148   //
    149   if (fSPhi+fDPhi>twopi) fSPhi-=twopi;
    150 
    151   // Check theta angles
    152 
    153   if (pSTheta<0 || pSTheta>pi)
    154   {
    155     G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
    156     G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
    157                 "stheta outside 0-PI range.");
    158   }
    159   else
    160   {
    161     fSTheta=pSTheta;
    162   }
    163 
    164   if (pDTheta+pSTheta>=pi)
    165   {
    166     fDTheta=pi-pSTheta;
    167   }
    168   else if (pDTheta>0)
    169   {
    170     fDTheta=pDTheta;
    171   }
    172   else
    173   {
    174     G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
    175            << "        Negative delta-Theta ! - "
    176            << pDTheta << G4endl;
    177     G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
    178                 "Invalid pDTheta.");
    179   }
     121  CheckPhiAngles(pSPhi, pDPhi);
     122  CheckThetaAngles(pSTheta, pDTheta);
    180123}
    181124
    … …  
    219162                                        G4double& pMin, G4double& pMax ) const
    220163{
    221   if ( fDPhi==twopi && fDTheta==pi)  // !pTransform.IsRotated() &&
     164  if ( fFullSphere )
    222165  {
    223166    // Special case handling for solid spheres-shells
    … …  
    310253        yoff1=yoffset-yMin;
    311254        yoff2=yMax-yoffset;
    312         if (yoff1>=0&&yoff2>=0)
     255        if ((yoff1>=0) && (yoff2>=0))
    313256        {
    314257          // Y limits cross max/min x => no change
    … …  
    334277        xoff1=xoffset-xMin;
    335278        xoff2=xMax-xoffset;
    336         if (xoff1>=0&&xoff2>=0)
     279        if ((xoff1>=0) && (xoff2>=0))
    337280        {
    338281          // X limits cross max/min y => no change
    … …  
    416359    }
    417360     
    418     if (pMin!=kInfinity || pMax!=-kInfinity)
     361    if ((pMin!=kInfinity) || (pMax!=-kInfinity))
    419362    {
    420363      existsAfterClip=true;
    … …  
    453396// Return whether point inside/outside/on surface
    454397// Split into radius, phi, theta checks
    455 // Each check modifies `in', or returns as approprate
     398// Each check modifies 'in', or returns as approprate
    456399
    457400EInside G4Sphere::Inside( const G4ThreeVector& p ) const
    … …  
    459402  G4double rho,rho2,rad2,tolRMin,tolRMax;
    460403  G4double pPhi,pTheta;
    461   EInside in=kOutside;
     404  EInside in = kOutside;
     405  static const G4double halfAngTolerance = kAngTolerance*0.5;
     406  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
     407  const G4double halfRminTolerance = fRminTolerance*0.5;
     408  const G4double Rmax_minus = fRmax - halfRmaxTolerance;
     409  const G4double Rmin_plus  = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
    462410
    463411  rho2 = p.x()*p.x() + p.y()*p.y() ;
    464412  rad2 = rho2 + p.z()*p.z() ;
    465413
    466   //  if(rad2 >= 1.369e+19) DBG();
    467   //  G4double rad = std::sqrt(rad2);
    468   // Check radial surfaces
    469   // sets `in'
    470 
    471   if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5;
    472   else         tolRMin = 0 ;
    473  
    474   tolRMax = fRmax - kRadTolerance*0.5 ;
    475   //  const G4double  fractionTolerance = 1.0e-12;
    476   const G4double  flexRadMaxTolerance = // kRadTolerance;
    477     std::max(kRadTolerance, fEpsilon * fRmax);
    478 
    479   const G4double  Rmax_minus = fRmax - flexRadMaxTolerance*0.5;
    480   const G4double  flexRadMinTolerance = std::max(kRadTolerance,
    481                      fEpsilon * fRmin);
    482   const G4double  Rmin_plus = (fRmin > 0) ? fRmin + flexRadMinTolerance*0.5 : 0 ;
     414  // Check radial surfaces. Sets 'in'
     415
     416  tolRMin = Rmin_plus;
     417  tolRMax = Rmax_minus;
     418
     419  if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
     420  {
     421    in = kInside;
     422  }
     423  else
     424  {
     425    tolRMax = fRmax + halfRmaxTolerance;                  // outside case
     426    tolRMin = std::max(fRmin-halfRminTolerance, 0.);      // outside case
     427    if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
     428    {
     429      in = kSurface;
     430    }
     431    else
     432    {
     433      return in = kOutside;
     434    }
     435  }
     436
     437  // Phi boundaries   : Do not check if it has no phi boundary!
     438
     439  if ( !fFullPhiSphere && rho2 )  // [fDPhi < twopi] and [p.x or p.y]
     440  {
     441    pPhi = std::atan2(p.y(),p.x()) ;
     442
     443    if      ( pPhi < fSPhi - halfAngTolerance  ) { pPhi += twopi; }
     444    else if ( pPhi > ePhi + halfAngTolerance )   { pPhi -= twopi; }
    483445   
    484 if(rad2 <= Rmax_minus*Rmax_minus && rad2 >= Rmin_plus*Rmin_plus) in = kInside ;
    485 
    486 // if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kInside ;
    487   // if ( rad <= tolRMax && rad >= tolRMin )  in = kInside ;
    488   else
    489   {
    490     tolRMax = fRmax + kRadTolerance*0.5 ;
    491     tolRMin = fRmin - kRadTolerance*0.5 ;
    492 
    493     if ( tolRMin < 0.0 ) tolRMin = 0.0 ;
    494    
    495      if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kSurface ;
    496     //  if ( rad <= tolRMax && rad >= tolRMin )  in = kSurface ;
    497     else                                                return in = kOutside ;
    498   }
    499 
    500   // Phi boundaries   : Do not check if it has no phi boundary!
    501   // (in != kOutside). It is new J.Apostolakis proposal of 30.10.03
    502 
    503   if ( ( fDPhi < twopi - kAngTolerance ) &&
    504        ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
    505   {
    506     pPhi = std::atan2(p.y(),p.x()) ;
    507 
    508     if      ( pPhi < fSPhi - kAngTolerance*0.5  )         pPhi += twopi ;
    509     else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 )  pPhi -= twopi;
    510    
    511     if ((pPhi < fSPhi - kAngTolerance*0.5) || 
    512         (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) )  return in = kOutside ;
     446    if ( (pPhi < fSPhi - halfAngTolerance)
     447      || (pPhi > ePhi + halfAngTolerance) )      { return in = kOutside; }
    513448   
    514449    else if (in == kInside)  // else it's kSurface anyway already
    515450    {
    516       if ( (pPhi < fSPhi + kAngTolerance*0.5) ||
    517            (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) )      in = kSurface ;       
     451      if ( (pPhi < fSPhi + halfAngTolerance)
     452        || (pPhi > ePhi - halfAngTolerance) )    { in = kSurface; }     
    518453    }
    519454  }
    520455
    521456  // Theta bondaries
    522   // (in!=kOutside)
    523457 
    524   if ( (rho2 || p.z()) && fDTheta < pi - kAngTolerance*0.5 )
     458  if ( (rho2 || p.z()) && (!fFullThetaSphere) )
    525459  {
    526460    rho    = std::sqrt(rho2);
    … …  
    529463    if ( in == kInside )
    530464    {
    531       if ( (pTheta < fSTheta + kAngTolerance*0.5)
    532         || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) )
    533       {
    534         if ( (pTheta >= fSTheta - kAngTolerance*0.5)
    535           && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) )
    536         {
    537           in = kSurface ;
     465      if ( (pTheta < fSTheta + halfAngTolerance)
     466        || (pTheta > eTheta - halfAngTolerance) )
     467      {
     468        if ( (pTheta >= fSTheta - halfAngTolerance)
     469          && (pTheta <= eTheta + halfAngTolerance) )
     470        {
     471          in = kSurface;
    538472        }
    539473        else
    540474        {
    541           in = kOutside ;
     475          in = kOutside;
    542476        }
    543477      }
    … …  
    545479    else
    546480    {
    547       if ( (pTheta < fSTheta - kAngTolerance*0.5)
    548         || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) )
    549       {
    550         in = kOutside ;
     481      if ( (pTheta < fSTheta - halfAngTolerance)
     482        || (pTheta > eTheta + halfAngTolerance) )
     483      {
     484        in = kOutside;
    551485      }
    552486    }
    … …  
    568502  G4double distSPhi = kInfinity, distEPhi = kInfinity;
    569503  G4double distSTheta = kInfinity, distETheta = kInfinity;
    570   G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
    571504  G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
    572505  G4ThreeVector norm, sumnorm(0.,0.,0.);
     506
     507  static const G4double halfCarTolerance = 0.5*kCarTolerance;
     508  static const G4double halfAngTolerance = 0.5*kAngTolerance;
    573509
    574510  rho2 = p.x()*p.x()+p.y()*p.y();
    … …  
    579515  if (fRmin)  distRMin = std::fabs(rad-fRmin);
    580516   
    581   if ( rho && (fDPhi < twopi || fDTheta < pi) )
     517  if ( rho && !fFullSphere )
    582518  {
    583519    pPhi = std::atan2(p.y(),p.x());
    584520
    585     if(pPhi  < fSPhi-dAngle)           pPhi     += twopi;
    586     else if(pPhi > fSPhi+fDPhi+dAngle) pPhi     -= twopi;
    587   }
    588   if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
     521    if (pPhi < fSPhi-halfAngTolerance)     { pPhi += twopi; }
     522    else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
     523  }
     524  if ( !fFullPhiSphere )
    589525  {
    590526    if ( rho )
    591527    {
    592       distSPhi = std::fabs( pPhi - fSPhi );
    593       distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
     528      distSPhi = std::fabs( pPhi-fSPhi );
     529      distEPhi = std::fabs( pPhi-ePhi );
    594530    }
    595531    else if( !fRmin )
    … …  
    598534      distEPhi = 0.;
    599535    }
    600     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
    601     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
     536    nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
     537    nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
    602538  }       
    603   if ( fDTheta < pi ) // && rad ) // old limitation against (0,0,0)
     539  if ( !fFullThetaSphere )
    604540  {
    605541    if ( rho )
    … …  
    607543      pTheta     = std::atan2(rho,p.z());
    608544      distSTheta = std::fabs(pTheta-fSTheta);
    609       distETheta = std::fabs(pTheta-fSTheta-fDTheta);
     545      distETheta = std::fabs(pTheta-eTheta);
    610546 
    611       nTs = G4ThreeVector(-std::cos(fSTheta)*p.x()/rho, //  *std::cos(pPhi),
    612                           -std::cos(fSTheta)*p.y()/rho, //  *std::sin(pPhi),
    613                            std::sin(fSTheta)                   );
    614 
    615       nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*p.x()/rho, // *std::cos(pPhi),
    616                            std::cos(fSTheta+fDTheta)*p.y()/rho, // *std::sin(pPhi),
    617                           -std::sin(fSTheta+fDTheta)                  );   
     547      nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
     548                          -cosSTheta*p.y()/rho,
     549                           sinSTheta          );
     550
     551      nTe = G4ThreeVector( cosETheta*p.x()/rho,
     552                           cosETheta*p.y()/rho,
     553                          -sinETheta          );   
    618554    }
    619555    else if( !fRmin )
    … …  
    622558      {             
    623559        distSTheta = 0.;
    624         nTs = G4ThreeVector(0.,0.,-1.);
    625       }
    626       if ( fSTheta + fDTheta < pi ) // distETheta = 0.;
     560        nTs = G4ThreeVector(0.,0.,-1.);
     561      }
     562      if ( eTheta < pi )
    627563      {             
    628564        distETheta = 0.;
    629         nTe = G4ThreeVector(0.,0.,1.);
     565        nTe = G4ThreeVector(0.,0.,1.);
    630566      }
    631567    }   
    632568  }
    633   if( rad )  nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
    634 
    635   if( distRMax <= delta )
     569  if( rad )  { nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad); }
     570
     571  if( distRMax <= halfCarTolerance )
    636572  {
    637573    noSurfaces ++;
    638574    sumnorm += nR;
    639575  }
    640   if( fRmin && distRMin <= delta )
     576  if( fRmin && (distRMin <= halfCarTolerance) )
    641577  {
    642578    noSurfaces ++;
    643579    sumnorm -= nR;
    644580  }
    645   if( fDPhi < twopi )   
    646   {
    647     if (distSPhi <= dAngle)
     581  if( !fFullPhiSphere )   
     582  {
     583    if (distSPhi <= halfAngTolerance)
    648584    {
    649585      noSurfaces ++;
    650586      sumnorm += nPs;
    651587    }
    652     if (distEPhi <= dAngle)
     588    if (distEPhi <= halfAngTolerance)
    653589    {
    654590      noSurfaces ++;
    … …  
    656592    }
    657593  }
    658   if ( fDTheta < pi )
    659   {
    660     if (distSTheta <= dAngle && fSTheta > 0.)
     594  if ( !fFullThetaSphere )
     595  {
     596    if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
    661597    {
    662598      noSurfaces ++;
    663       if( rad <= delta && fDPhi >= twopi) sumnorm += nZ;
    664       else                                sumnorm += nTs;
    665     }
    666     if (distETheta <= dAngle && fSTheta+fDTheta < pi)
     599      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm += nZ;  }
     600      else                                              { sumnorm += nTs; }
     601    }
     602    if ((distETheta <= halfAngTolerance) && (eTheta < pi))
    667603    {
    668604      noSurfaces ++;
    669       if( rad <= delta && fDPhi >= twopi) sumnorm -= nZ;
    670       else                                sumnorm += nTe;
    671       if(sumnorm.z() == 0.)               sumnorm += nZ;
     605      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm -= nZ;  }
     606      else                                              { sumnorm += nTe; }
     607      if(sumnorm.z() == 0.)  { sumnorm += nZ; }
    672608    }
    673609  }
    … …  
    680616     norm = ApproxSurfaceNormal(p);
    681617  }
    682   else if ( noSurfaces == 1 ) norm = sumnorm;
    683   else                        norm = sumnorm.unit();
     618  else if ( noSurfaces == 1 )  { norm = sumnorm; }
     619  else                         { norm = sumnorm.unit(); }
    684620  return norm;
    685621}
    … …  
    735671   
    736672  pPhi = std::atan2(p.y(),p.x());
    737   if (pPhi<0) pPhi += twopi;
    738 
    739   if (fDPhi<twopi&&rho)
     673  if (pPhi<0) { pPhi += twopi; }
     674
     675  if (!fFullPhiSphere && rho)
    740676  {
    741677    if (fSPhi<0)
    … …  
    774710  //
    775711
    776   if (fDTheta<pi&&rad)
     712  if (!fFullThetaSphere && rad)
    777713  {
    778714    pTheta=std::atan2(rho,p.z());
    … …  
    809745      break;
    810746    case kNSPhi:
    811       norm=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
     747      norm=G4ThreeVector(sinSPhi,-cosSPhi,0);
    812748      break;
    813749    case kNEPhi:
    814       norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
     750      norm=G4ThreeVector(-sinEPhi,cosEPhi,0);
    815751      break;
    816752    case kNSTheta:
    817       norm=G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
    818                          -std::cos(fSTheta)*std::sin(pPhi),
    819                           std::sin(fSTheta)            );
    820       //  G4cout<<G4endl<<" case kNSTheta:"<<G4endl;
    821       //  G4cout<<"pPhi = "<<pPhi<<G4endl;
    822       //  G4cout<<"rad  = "<<rad<<G4endl;
    823       //  G4cout<<"pho  = "<<rho<<G4endl;
    824       //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
    825       //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
     753      norm=G4ThreeVector(-cosSTheta*std::cos(pPhi),
     754                         -cosSTheta*std::sin(pPhi),
     755                          sinSTheta            );
    826756      break;
    827757    case kNETheta:
    828       norm=G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
    829                           std::cos(fSTheta+fDTheta)*std::sin(pPhi),
    830                          -std::sin(fSTheta+fDTheta)              );
    831 
    832       //  G4cout<<G4endl<<" case kNETheta:"<<G4endl;
    833       //  G4cout<<"pPhi = "<<pPhi<<G4endl;
    834       //  G4cout<<"rad  = "<<rad<<G4endl;
    835       //  G4cout<<"pho  = "<<rho<<G4endl;
    836       //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
    837       //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
     758      norm=G4ThreeVector( cosETheta*std::cos(pPhi),
     759                          cosETheta*std::sin(pPhi),
     760                         -sinETheta              );
    838761      break;
    839762    default:
    840763      DumpInfo();
    841       G4Exception("G4Sphere::ApproxSurfaceNormal()", "Notification", JustWarning,
     764      G4Exception("G4Sphere::ApproxSurfaceNormal()","Notification",JustWarning,
    842765                  "Undefined side for valid surface normal to solid.");
    843766      break;   
    844   } // end case
     767  }
    845768
    846769  return norm;
    … …  
    880803{
    881804  G4double snxt = kInfinity ;      // snxt = default return value
    882 
    883805  G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
    884 
    885   G4double tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ;
    886806  G4double tolSTheta=0., tolETheta=0. ;
     807  const G4double dRmax = 100.*fRmax;
     808
     809  static const G4double halfCarTolerance = kCarTolerance*0.5;
     810  static const G4double halfAngTolerance = kAngTolerance*0.5;
     811  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
     812  const G4double halfRminTolerance = fRminTolerance*0.5;
     813  const G4double tolORMin2 = (fRmin>halfRminTolerance)
     814               ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
     815  const G4double tolIRMin2 =
     816               (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
     817  const G4double tolORMax2 =
     818               (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
     819  const G4double tolIRMax2 =
     820               (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
    887821
    888822  // Intersection point
    889 
     823  //
    890824  G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
    891825
    892826  // Phi intersection
    893 
    894   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ;
    895 
    896   // Phi flag and precalcs
    897 
    898   G4bool segPhi ;       
    899   G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ;
    900   G4double cosHDPhiOT=0., cosHDPhiIT=0. ;
     827  //
     828  G4double Comp ;
     829
     830  // Phi precalcs
     831  //
    901832  G4double Dist, cosPsi ;
    902833
    903   G4bool segTheta ;                             // Theta flag and precals
    904   G4double tanSTheta, tanETheta ;
    905   G4double tanSTheta2, tanETheta2 ;
     834  // Theta precalcs
     835  //
    906836  G4double dist2STheta, dist2ETheta ;
    907837  G4double t1, t2, b, c, d2, d, s = kInfinity ;
    908838
    909839  // General Precalcs
    910 
     840  //
    911841  rho2 = p.x()*p.x() + p.y()*p.y() ;
    912842  rad2 = rho2 + p.z()*p.z() ;
    … …  
    916846  pDotV3d = pDotV2d + p.z()*v.z() ;
    917847
    918   // Radial Precalcs
    919 
    920   if (fRmin > kRadTolerance*0.5)
    921   {
    922     tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5);
    923   }
    924   else
    925   {
    926     tolORMin2 = 0 ;
    927   }
    928   tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ;
    929   tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ;
    930   tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ;
    931 
    932   // Set phi divided flag and precalcs
    933 
    934   if (fDPhi < twopi)
    935   {
    936     segPhi = true ;
    937     hDPhi = 0.5*fDPhi ;    // half delta phi
    938     cPhi = fSPhi + hDPhi ;
    939 
    940     hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
    941     hDPhiIT = hDPhi-0.5*kAngTolerance;
    942 
    943     sinCPhi    = std::sin(cPhi) ;
    944     cosCPhi    = std::cos(cPhi) ;
    945     cosHDPhiOT = std::cos(hDPhiOT) ;
    946     cosHDPhiIT = std::cos(hDPhiIT) ;
    947   }
    948   else
    949   {
    950     segPhi = false ;
    951   }
    952 
    953848  // Theta precalcs
    954    
    955   if (fDTheta < pi )
    956   {
    957     segTheta  = true ;
    958     tolSTheta = fSTheta - kAngTolerance*0.5 ;
    959     tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ;
    960   }
    961   else
    962   {
    963     segTheta = false ;
     849  //
     850  if (!fFullThetaSphere)
     851  {
     852    tolSTheta = fSTheta - halfAngTolerance ;
     853    tolETheta = eTheta + halfAngTolerance ;
    964854  }
    965855
    … …  
    979869
    980870  c = rad2 - fRmax*fRmax ;
    981   const G4double  flexRadMaxTolerance = // kRadTolerance;
    982     std::max(kRadTolerance, fEpsilon * fRmax);
    983 
    984   //  if (c > kRadTolerance*fRmax)
    985   if (c > flexRadMaxTolerance*fRmax)
    986   {
    987     // If outside toleranct boundary of outer G4Sphere
    988     // [should be std::sqrt(rad2)-fRmax > kRadTolerance*0.5]
     871
     872  if (c > fRmaxTolerance*fRmax)
     873  {
     874    // If outside tolerant boundary of outer G4Sphere
     875    // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
    989876
    990877    d2 = pDotV3d*pDotV3d - c ;
    … …  
    996883      if (s >= 0 )
    997884      {
     885        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
     886        {              // 64 bits systems. Split long distances and recompute
     887          G4double fTerm = s-std::fmod(s,dRmax);
     888          s = fTerm + DistanceToIn(p+fTerm*v,v);
     889        }
    998890        xi   = p.x() + s*v.x() ;
    999891        yi   = p.y() + s*v.y() ;
    1000892        rhoi = std::sqrt(xi*xi + yi*yi) ;
    1001893
    1002         if (segPhi && rhoi)    // Check phi intersection
     894        if (!fFullPhiSphere && rhoi)    // Check phi intersection
    1003895        {
    1004896          cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
    … …  
    1006898          if (cosPsi >= cosHDPhiOT)
    1007899          {
    1008             if (segTheta)   // Check theta intersection
     900            if (!fFullThetaSphere)   // Check theta intersection
    1009901            {
    1010902              zi = p.z() + s*v.z() ;
    … …  
    1027919        else
    1028920        {
    1029           if (segTheta)    // Check theta intersection
     921          if (!fFullThetaSphere)    // Check theta intersection
    1030922          {
    1031923            zi = p.z() + s*v.z() ;
    … …  
    1059951    d2 = pDotV3d*pDotV3d - c ;
    1060952
    1061     // if (rad2 > tolIRMin2 && pDotV3d < 0 )
    1062 
    1063     if (rad2 > tolIRMax2 && ( d2 >= flexRadMaxTolerance*fRmax && pDotV3d < 0 ) )
    1064     {
    1065       if (segPhi)
     953    if ( (rad2 > tolIRMax2)
     954      && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
     955    {
     956      if (!fFullPhiSphere)
    1066957      {
    1067958        // Use inner phi tolerant boundary -> if on tolerant
    … …  
    1074965          // inside radii, delta r -ve, inside phi
    1075966
    1076           if (segTheta)
     967          if ( !fFullThetaSphere )
    1077968          {
    1078969            if ( (pTheta >= tolSTheta + kAngTolerance)
    … …  
    1090981      else
    1091982      {
    1092         if ( segTheta )
     983        if ( !fFullThetaSphere )
    1093984        {
    1094985          if ( (pTheta >= tolSTheta + kAngTolerance)
    … …  
    11091000  // - Always farthest root, because would have passed through outer
    11101001  //   surface first.
    1111   // - Tolerant check for if travelling through solid
     1002  // - Tolerant check if travelling through solid
    11121003
    11131004  if (fRmin)
    … …  
    11191010    // Check for immediate entry/already inside and travelling outwards
    11201011
    1121     // if (c >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMin2 )
    1122 
    1123     if ( c > -kRadTolerance*0.5 && rad2 < tolIRMin2 &&
    1124          ( d2 < fRmin*kCarTolerance || pDotV3d >= 0 ) )
    1125     {
    1126       if (segPhi)
     1012    if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
     1013      && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
     1014    {
     1015      if ( !fFullPhiSphere )
    11271016      {
    11281017        // Use inner phi tolerant boundary -> if on tolerant
    … …  
    11341023          // inside radii, delta r -ve, inside phi
    11351024          //
    1136           if (segTheta)
     1025          if ( !fFullThetaSphere )
    11371026          {
    11381027            if ( (pTheta >= tolSTheta + kAngTolerance)
    … …  
    11501039      else
    11511040      {
    1152         if (segTheta)
     1041        if ( !fFullThetaSphere )
    11531042        {
    11541043          if ( (pTheta >= tolSTheta + kAngTolerance)
    … …  
    11661055    else   // Not special tolerant case
    11671056    {
    1168       //  d2 = pDotV3d*pDotV3d - c ;
    1169 
    11701057      if (d2 >= 0)
    11711058      {
    11721059        s = -pDotV3d + std::sqrt(d2) ;
    1173         if ( s >= kRadTolerance*0.5 )  // It was >= 0 ??
     1060        if ( s >= halfRminTolerance )  // It was >= 0 ??
    11741061        {
    11751062          xi   = p.x() + s*v.x() ;
    … …  
    11771064          rhoi = std::sqrt(xi*xi+yi*yi) ;
    11781065
    1179           if ( segPhi && rhoi )   // Check phi intersection
     1066          if ( !fFullPhiSphere && rhoi )   // Check phi intersection
    11801067          {
    11811068            cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
    … …  
    11831070            if (cosPsi >= cosHDPhiOT)
    11841071            {
    1185               if (segTheta)  // Check theta intersection
     1072              if ( !fFullThetaSphere )  // Check theta intersection
    11861073              {
    11871074                zi = p.z() + s*v.z() ;
    … …  
    12041091          else
    12051092          {
    1206             if (segTheta)   // Check theta intersection
     1093            if ( !fFullThetaSphere )   // Check theta intersection
    12071094            {
    12081095              zi = p.z() + s*v.z() ;
    … …  
    12141101              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
    12151102              {
    1216                 snxt = s ;
     1103                snxt = s;
    12171104              }
    12181105            }
    12191106            else
    12201107            {
    1221               snxt=s;
     1108              snxt = s;
    12221109            }
    12231110          }
    … …  
    12361123  //         -> Should use some form of loop Construct
    12371124  //
    1238   if ( segPhi )
    1239   {
    1240     // First phi surface (`S'tarting phi)
    1241 
    1242     sinSPhi = std::sin(fSPhi) ;
    1243     cosSPhi = std::cos(fSPhi) ;
    1244 
     1125  if ( !fFullPhiSphere )
     1126  {
     1127    // First phi surface ('S'tarting phi)
    12451128    // Comp = Component in outwards normal dirn
    12461129    //
    1247     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
     1130    Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
    12481131                   
    12491132    if ( Comp < 0 )
    … …  
    12511134      Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
    12521135
    1253       if (Dist < kCarTolerance*0.5)
     1136      if (Dist < halfCarTolerance)
    12541137      {
    12551138        s = Dist/Comp ;
    … …  
    12821165            // (=>intersect at origin =>fRmax=0)
    12831166            //
    1284             if ( segTheta )
     1167            if ( !fFullThetaSphere )
    12851168            {
    12861169              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
    … …  
    13051188    }
    13061189
    1307     // Second phi surface (`E'nding phi)
    1308 
    1309     ePhi    = fSPhi + fDPhi ;
    1310     sinEPhi = std::sin(ePhi)     ;
    1311     cosEPhi = std::cos(ePhi)     ;
    1312 
    1313     // Compnent in outwards normal dirn
    1314 
    1315     Comp    = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
     1190    // Second phi surface ('E'nding phi)
     1191    // Component in outwards normal dirn
     1192
     1193    Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
    13161194       
    13171195    if (Comp < 0)
    13181196    {
    13191197      Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
    1320       if ( Dist < kCarTolerance*0.5 )
     1198      if ( Dist < halfCarTolerance )
    13211199      {
    13221200        s = Dist/Comp ;
    … …  
    13401218            rhoi2 = rho2  ;
    13411219            radi2 = rad2  ;
    1342           } if ( (radi2 <= tolORMax2)
     1220          }
     1221          if ( (radi2 <= tolORMax2)
    13431222            && (radi2 >= tolORMin2)
    13441223            && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
    … …  
    13481227            // (=>intersect at origin =>fRmax=0)
    13491228            //
    1350             if ( segTheta )
     1229            if ( !fFullThetaSphere )
    13511230            {
    13521231              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
    … …  
    13741253  // Theta segment intersection
    13751254
    1376   if ( segTheta )
     1255  if ( !fFullThetaSphere )
    13771256  {
    13781257
    … …  
    13961275    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
    13971276
    1398     tanSTheta  = std::tan(fSTheta)         ;
    1399     tanSTheta2 = tanSTheta*tanSTheta  ;
    1400     tanETheta  = std::tan(fSTheta+fDTheta) ;
    1401     tanETheta2 = tanETheta*tanETheta  ;
    1402      
    14031277    if (fSTheta)
    14041278    {
    … …  
    14091283      dist2STheta = kInfinity ;
    14101284    }
    1411     if ( fSTheta + fDTheta < pi )
     1285    if ( eTheta < pi )
    14121286    {
    14131287      dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
    14141288    }
    1415       else
     1289    else
    14161290    {
    14171291      dist2ETheta=kInfinity;
    14181292    }     
    1419     if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0)
     1293    if ( pTheta < tolSTheta )
    14201294    {
    14211295      // Inside (theta<stheta-tol) s theta cone
    … …  
    14241298      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
    14251299      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
    1426        
    1427       b  = t2/t1 ;
    1428       c  = dist2STheta/t1 ;
    1429       d2 = b*b - c ;
    1430 
    1431       if ( d2 >= 0 )
    1432       {
    1433         d = std::sqrt(d2) ;
    1434         s = -b - d ;    // First root
    1435         zi    = p.z() + s*v.z();
    1436 
    1437         if ( s < 0 || zi*(fSTheta - halfpi) > 0 )
    1438         {
    1439           s = -b+d;    // Second root
    1440         }
    1441         if (s >= 0 && s < snxt)
    1442         {
    1443           xi    = p.x() + s*v.x();
    1444           yi    = p.y() + s*v.y();
     1300      if (t1)
     1301      {   
     1302        b  = t2/t1 ;
     1303        c  = dist2STheta/t1 ;
     1304        d2 = b*b - c ;
     1305
     1306        if ( d2 >= 0 )
     1307        {
     1308          d = std::sqrt(d2) ;
     1309          s = -b - d ;    // First root
    14451310          zi    = p.z() + s*v.z();
    1446           rhoi2 = xi*xi + yi*yi;
    1447           radi2 = rhoi2 + zi*zi;
    1448           if ( (radi2 <= tolORMax2)
    1449             && (radi2 >= tolORMin2)
    1450             && (zi*(fSTheta - halfpi) <= 0) )
    1451           {
    1452             if ( segPhi && rhoi2 )  // Check phi intersection
    1453             {
    1454               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1455               if (cosPsi >= cosHDPhiOT)
     1311
     1312          if ( (s < 0) || (zi*(fSTheta - halfpi) > 0) )
     1313          {
     1314            s = -b+d;    // Second root
     1315          }
     1316          if ((s >= 0) && (s < snxt))
     1317          {
     1318            xi    = p.x() + s*v.x();
     1319            yi    = p.y() + s*v.y();
     1320            zi    = p.z() + s*v.z();
     1321            rhoi2 = xi*xi + yi*yi;
     1322            radi2 = rhoi2 + zi*zi;
     1323            if ( (radi2 <= tolORMax2)
     1324              && (radi2 >= tolORMin2)
     1325              && (zi*(fSTheta - halfpi) <= 0) )
     1326            {
     1327              if ( !fFullPhiSphere && rhoi2 )  // Check phi intersection
     1328              {
     1329                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1330                if (cosPsi >= cosHDPhiOT)
     1331                {
     1332                  snxt = s ;
     1333                }
     1334              }
     1335              else
    14561336              {
    14571337                snxt = s ;
    14581338              }
    1459             }
    1460             else
    1461             {
    1462               snxt = s ;
    14631339            }
    14641340          }
    … …  
    14691345      // Second >= 0 root should be considered
    14701346       
    1471       if ( fSTheta + fDTheta < pi )
     1347      if ( eTheta < pi )
    14721348      {
    14731349        t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
    14741350        t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
    1475        
     1351        if (t1)
     1352        {
     1353          b  = t2/t1 ;
     1354          c  = dist2ETheta/t1 ;
     1355          d2 = b*b - c ;
     1356
     1357          if (d2 >= 0)
     1358          {
     1359            d = std::sqrt(d2) ;
     1360            s = -b + d ;    // Second root
     1361
     1362            if ( (s >= 0) && (s < snxt) )
     1363            {
     1364              xi    = p.x() + s*v.x() ;
     1365              yi    = p.y() + s*v.y() ;
     1366              zi    = p.z() + s*v.z() ;
     1367              rhoi2 = xi*xi + yi*yi   ;
     1368              radi2 = rhoi2 + zi*zi   ;
     1369
     1370              if ( (radi2 <= tolORMax2)
     1371                && (radi2 >= tolORMin2)
     1372                && (zi*(eTheta - halfpi) <= 0) )
     1373              {
     1374                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
     1375                {
     1376                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1377                  if (cosPsi >= cosHDPhiOT)
     1378                  {
     1379                    snxt = s ;
     1380                  }
     1381                }
     1382                else
     1383                {
     1384                  snxt = s ;
     1385                }
     1386              }
     1387            }
     1388          }
     1389        }
     1390      }
     1391    } 
     1392    else if ( pTheta > tolETheta )
     1393    {
     1394      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
     1395      // Inside (theta > etheta+tol) e-theta cone
     1396      // First root of etheta cone, second if first root 'imaginary'
     1397
     1398      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
     1399      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
     1400      if (t1)
     1401      { 
    14761402        b  = t2/t1 ;
    14771403        c  = dist2ETheta/t1 ;
    … …  
    14811407        {
    14821408          d = std::sqrt(d2) ;
    1483           s = -b + d ;    // Second root
    1484 
    1485           if (s >= 0 && s < snxt)
     1409          s = -b - d ;    // First root
     1410          zi    = p.z() + s*v.z();
     1411
     1412          if ( (s < 0) || (zi*(eTheta - halfpi) > 0) )
     1413          {
     1414            s = -b + d ;           // second root
     1415          }
     1416          if ( (s >= 0) && (s < snxt) )
    14861417          {
    14871418            xi    = p.x() + s*v.x() ;
    … …  
    14921423
    14931424            if ( (radi2 <= tolORMax2)
    1494               && (radi2 >= tolORMin2)
    1495               && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
    1496             {
    1497               if (segPhi && rhoi2)   // Check phi intersection
     1425              && (radi2 >= tolORMin2) 
     1426              && (zi*(eTheta - halfpi) <= 0) )
     1427            {
     1428              if (!fFullPhiSphere && rhoi2)  // Check phi intersection
    14981429              {
    14991430                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    … …  
    15111442        }
    15121443      }
    1513     } 
    1514     else if ( pTheta > tolETheta )
    1515     {
    1516       // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
    1517       // Inside (theta > etheta+tol) e-theta cone
    1518       // First root of etheta cone, second if first root `imaginary'
    1519 
    1520       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
    1521       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
    1522        
    1523       b  = t2/t1 ;
    1524       c  = dist2ETheta/t1 ;
    1525       d2 = b*b - c ;
    1526 
    1527       if (d2 >= 0)
    1528       {
    1529         d = std::sqrt(d2) ;
    1530         s = -b - d ;    // First root
    1531         zi    = p.z() + s*v.z();
    1532 
    1533         if (s < 0 || zi*(fSTheta + fDTheta - halfpi) > 0)
    1534         {
    1535           s = -b + d ;           // second root
    1536         }
    1537         if (s >= 0 && s < snxt)
    1538         {
    1539           xi    = p.x() + s*v.x() ;
    1540           yi    = p.y() + s*v.y() ;
    1541           zi    = p.z() + s*v.z() ;
    1542           rhoi2 = xi*xi + yi*yi   ;
    1543           radi2 = rhoi2 + zi*zi   ;
    1544 
    1545           if ( (radi2 <= tolORMax2)
    1546             && (radi2 >= tolORMin2)
    1547             && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
    1548           {
    1549             if (segPhi && rhoi2)  // Check phi intersection
    1550             {
    1551               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1552               if (cosPsi >= cosHDPhiOT)
    1553               {
    1554                 snxt = s ;
    1555               }
    1556             }
    1557             else
    1558             {
    1559               snxt = s ;
    1560             }
    1561           }
    1562         }
    1563       }
    15641444
    15651445      // Possible intersection with STheta cone.
    … …  
    15701450        t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
    15711451        t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
    1572 
     1452        if (t1)
     1453        {
     1454          b  = t2/t1 ;
     1455          c  = dist2STheta/t1 ;
     1456          d2 = b*b - c ;
     1457
     1458          if (d2 >= 0)
     1459          {
     1460            d = std::sqrt(d2) ;
     1461            s = -b + d ;    // Second root
     1462
     1463            if ( (s >= 0) && (s < snxt) )
     1464            {
     1465              xi    = p.x() + s*v.x() ;
     1466              yi    = p.y() + s*v.y() ;
     1467              zi    = p.z() + s*v.z() ;
     1468              rhoi2 = xi*xi + yi*yi   ;
     1469              radi2 = rhoi2 + zi*zi   ;
     1470
     1471              if ( (radi2 <= tolORMax2)
     1472                && (radi2 >= tolORMin2)
     1473                && (zi*(fSTheta - halfpi) <= 0) )
     1474              {
     1475                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
     1476                {
     1477                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1478                  if (cosPsi >= cosHDPhiOT)
     1479                  {
     1480                    snxt = s ;
     1481                  }
     1482                }
     1483                else
     1484                {
     1485                  snxt = s ;
     1486                }
     1487              }
     1488            }
     1489          }
     1490        }
     1491      } 
     1492    }     
     1493    else if ( (pTheta < tolSTheta + kAngTolerance)
     1494           && (fSTheta > halfAngTolerance) )
     1495    {
     1496      // In tolerance of stheta
     1497      // If entering through solid [r,phi] => 0 to in
     1498      // else try 2nd root
     1499
     1500      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
     1501      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
     1502        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
     1503        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
     1504      {
     1505        if (!fFullPhiSphere && rho2)  // Check phi intersection
     1506        {
     1507          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
     1508          if (cosPsi >= cosHDPhiIT)
     1509          {
     1510            return 0 ;
     1511          }
     1512        }
     1513        else
     1514        {
     1515          return 0 ;
     1516        }
     1517      }
     1518
     1519      // Not entering immediately/travelling through
     1520
     1521      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
     1522      if (t1)
     1523      {
    15731524        b  = t2/t1 ;
    15741525        c  = dist2STheta/t1 ;
    … …  
    15781529        {
    15791530          d = std::sqrt(d2) ;
    1580           s = -b + d ;    // Second root
    1581 
    1582           if ( (s >= 0) && (s < snxt) )
    1583           {
     1531          s = -b + d ;
     1532          if ( (s >= halfCarTolerance) && (s < snxt) && (fSTheta < halfpi) )
     1533          {  // ^^^^^^^^^^^^^^^^^^^^^  shouldn't it be >=0 instead ?
    15841534            xi    = p.x() + s*v.x() ;
    15851535            yi    = p.y() + s*v.y() ;
    … …  
    15921542              && (zi*(fSTheta - halfpi) <= 0) )
    15931543            {
    1594               if (segPhi && rhoi2)   // Check phi intersection
     1544              if ( !fFullPhiSphere && rhoi2 )    // Check phi intersection
     1545              {
     1546                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1547                if ( cosPsi >= cosHDPhiOT )
     1548                {
     1549                  snxt = s ;
     1550                }
     1551              }
     1552              else
     1553              {
     1554                snxt = s ;
     1555              }
     1556            }
     1557          }
     1558        }
     1559      }
     1560    }   
     1561    else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
     1562    {
     1563
     1564      // In tolerance of etheta
     1565      // If entering through solid [r,phi] => 0 to in
     1566      // else try 2nd root
     1567
     1568      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
     1569
     1570      if (   ((t2<0) && (eTheta < halfpi)
     1571          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
     1572        ||   ((t2>=0) && (eTheta > halfpi)
     1573          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
     1574        ||   ((v.z()>0) && (eTheta == halfpi)
     1575          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))  )
     1576      {
     1577        if (!fFullPhiSphere && rho2)   // Check phi intersection
     1578        {
     1579          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
     1580          if (cosPsi >= cosHDPhiIT)
     1581          {
     1582            return 0 ;
     1583          }
     1584        }
     1585        else
     1586        {
     1587          return 0 ;
     1588        }
     1589      }
     1590
     1591      // Not entering immediately/travelling through
     1592
     1593      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
     1594      if (t1)
     1595      {
     1596        b  = t2/t1 ;
     1597        c  = dist2ETheta/t1 ;
     1598        d2 = b*b - c ;
     1599
     1600        if (d2 >= 0)
     1601        {
     1602          d = std::sqrt(d2) ;
     1603          s = -b + d ;
     1604       
     1605          if ( (s >= halfCarTolerance)
     1606            && (s < snxt) && (eTheta > halfpi) )
     1607          {
     1608            xi    = p.x() + s*v.x() ;
     1609            yi    = p.y() + s*v.y() ;
     1610            zi    = p.z() + s*v.z() ;
     1611            rhoi2 = xi*xi + yi*yi   ;
     1612            radi2 = rhoi2 + zi*zi   ;
     1613
     1614            if ( (radi2 <= tolORMax2)
     1615              && (radi2 >= tolORMin2)
     1616              && (zi*(eTheta - halfpi) <= 0) )
     1617            {
     1618              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
    15951619              {
    15961620                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    … …  
    16061630            }
    16071631          }
    1608         }
    1609       } 
    1610     }     
    1611     else if ( (pTheta <tolSTheta + kAngTolerance)
    1612            && (fSTheta > kAngTolerance) )
    1613     {
    1614       // In tolerance of stheta
    1615       // If entering through solid [r,phi] => 0 to in
    1616       // else try 2nd root
    1617 
    1618       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
    1619       if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<pi*.5)
    1620         || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>pi*.5)
    1621         || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==pi*.5) )
    1622       {
    1623         if (segPhi && rho2)  // Check phi intersection
    1624         {
    1625           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
    1626           if (cosPsi >= cosHDPhiIT)
    1627           {
    1628             return 0 ;
    1629           }
    1630         }
    1631         else
    1632         {
    1633           return 0 ;
    1634         }
    1635       }
    1636 
    1637       // Not entering immediately/travelling through
    1638 
    1639       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
    1640       b  = t2/t1 ;
    1641       c  = dist2STheta/t1 ;
    1642       d2 = b*b - c ;
    1643 
    1644       if (d2 >= 0)
    1645       {
    1646         d = std::sqrt(d2) ;
    1647         s = -b + d ;
    1648         if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < pi*0.5) )
    1649         {
    1650           xi    = p.x() + s*v.x() ;
    1651           yi    = p.y() + s*v.y() ;
    1652           zi    = p.z() + s*v.z() ;
    1653           rhoi2 = xi*xi + yi*yi   ;
    1654           radi2 = rhoi2 + zi*zi   ;
    1655 
    1656           if ( (radi2 <= tolORMax2)
    1657             && (radi2 >= tolORMin2)
    1658             && (zi*(fSTheta - halfpi) <= 0) )
    1659           {
    1660             if ( segPhi && rhoi2 )    // Check phi intersection
    1661             {
    1662               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1663               if ( cosPsi >= cosHDPhiOT )
    1664               {
    1665                 snxt = s ;
    1666               }
    1667             }
    1668             else
    1669             {
    1670               snxt = s ;
    1671             }
    1672           }
    1673         }
    1674       }
    1675     }   
    1676     else if ( (pTheta > tolETheta - kAngTolerance)
    1677            && ((fSTheta + fDTheta) < pi-kAngTolerance) )   
    1678     {
    1679 
    1680       // In tolerance of etheta
    1681       // If entering through solid [r,phi] => 0 to in
    1682       // else try 2nd root
    1683 
    1684       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
    1685 
    1686       if (
    1687     (t2<0    && (fSTheta+fDTheta) <pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
    1688  || (t2>=0   && (fSTheta+fDTheta) >pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
    1689  || (v.z()>0 && (fSTheta+fDTheta)==pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
    1690          )
    1691       {
    1692         if (segPhi && rho2)   // Check phi intersection
    1693         {
    1694           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
    1695           if (cosPsi >= cosHDPhiIT)
    1696           {
    1697             return 0 ;
    1698           }
    1699         }
    1700         else
    1701         {
    1702           return 0 ;
    1703         }
    1704       }
    1705 
    1706       // Not entering immediately/travelling through
    1707 
    1708       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
    1709       b  = t2/t1 ;
    1710       c  = dist2ETheta/t1 ;
    1711       d2 = b*b - c ;
    1712 
    1713       if (d2 >= 0)
    1714       {
    1715         d = std::sqrt(d2) ;
    1716         s = -b + d ;
    1717        
    1718         if ( (s >= kCarTolerance*0.5)
    1719           && (s < snxt) && ((fSTheta + fDTheta) > pi*0.5) )
    1720         {
    1721           xi    = p.x() + s*v.x() ;
    1722           yi    = p.y() + s*v.y() ;
    1723           zi    = p.z() + s*v.z() ;
    1724           rhoi2 = xi*xi + yi*yi   ;
    1725           radi2 = rhoi2 + zi*zi   ;
    1726 
    1727           if ( (radi2 <= tolORMax2)
    1728             && (radi2 >= tolORMin2)
    1729             && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
    1730           {
    1731             if (segPhi && rhoi2)   // Check phi intersection
    1732             {
    1733               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1734               if (cosPsi>=cosHDPhiOT)
    1735               {
    1736                 snxt = s ;
    1737               }
    1738             }
    1739             else
    1740             {
    1741               snxt = s ;
    1742             }
    1743           }
    1744         }
    1745       }   
     1632        }
     1633      }   
    17461634    } 
    17471635    else
    … …  
    17521640      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
    17531641      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
    1754 
    1755       b  = t2/t1;
    1756       c  = dist2STheta/t1 ;
    1757       d2 = b*b - c ;
    1758 
    1759       if (d2 >= 0)
    1760       {
    1761         d = std::sqrt(d2) ;
    1762         s = -b + d ;    // second root
    1763 
    1764         if (s >= 0 && s < snxt)
    1765         {
    1766           xi    = p.x() + s*v.x() ;
    1767           yi    = p.y() + s*v.y() ;
    1768           zi    = p.z() + s*v.z() ;
    1769           rhoi2 = xi*xi + yi*yi   ;
    1770           radi2 = rhoi2 + zi*zi   ;
    1771 
    1772           if ( (radi2 <= tolORMax2)
    1773             && (radi2 >= tolORMin2)
    1774             && (zi*(fSTheta - halfpi) <= 0) )
    1775           {
    1776             if (segPhi && rhoi2)   // Check phi intersection
    1777             {
    1778               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1779               if (cosPsi >= cosHDPhiOT)
     1642      if (t1)
     1643      {
     1644        b  = t2/t1;
     1645        c  = dist2STheta/t1 ;
     1646        d2 = b*b - c ;
     1647
     1648        if (d2 >= 0)
     1649        {
     1650          d = std::sqrt(d2) ;
     1651          s = -b + d ;    // second root
     1652
     1653          if ((s >= 0) && (s < snxt))
     1654          {
     1655            xi    = p.x() + s*v.x() ;
     1656            yi    = p.y() + s*v.y() ;
     1657            zi    = p.z() + s*v.z() ;
     1658            rhoi2 = xi*xi + yi*yi   ;
     1659            radi2 = rhoi2 + zi*zi   ;
     1660
     1661            if ( (radi2 <= tolORMax2)
     1662              && (radi2 >= tolORMin2)
     1663              && (zi*(fSTheta - halfpi) <= 0) )
     1664            {
     1665              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
     1666              {
     1667                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1668                if (cosPsi >= cosHDPhiOT)
     1669                {
     1670                  snxt = s ;
     1671                }
     1672              }
     1673              else
    17801674              {
    17811675                snxt = s ;
    17821676              }
    1783             }
    1784             else
    1785             {
    1786               snxt = s ;
    17871677            }
    17881678          }
    … …  
    17911681      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
    17921682      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
    1793        
    1794       b  = t2/t1 ;
    1795       c  = dist2ETheta/t1 ;
    1796       d2 = b*b - c ;
    1797 
    1798       if (d2 >= 0)
    1799       {
    1800         d = std::sqrt(d2) ;
    1801         s = -b + d;    // second root
    1802 
    1803         if (s >= 0 && s < snxt)
    1804         {
    1805           xi    = p.x() + s*v.x() ;
    1806           yi    = p.y() + s*v.y() ;
    1807           zi    = p.z() + s*v.z() ;
    1808           rhoi2 = xi*xi + yi*yi   ;
    1809           radi2 = rhoi2 + zi*zi   ;
    1810 
    1811           if ( (radi2 <= tolORMax2)
    1812             && (radi2 >= tolORMin2)
    1813             && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
    1814           {
    1815             if (segPhi && rhoi2)   // Check phi intersection
    1816             {
    1817               cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
    1818               if ( cosPsi >= cosHDPhiOT )
    1819               {
    1820                 snxt=s;
    1821               }
    1822             }
    1823             else
    1824             {
    1825               snxt = s ;
     1683      if (t1)
     1684      {   
     1685        b  = t2/t1 ;
     1686        c  = dist2ETheta/t1 ;
     1687        d2 = b*b - c ;
     1688
     1689        if (d2 >= 0)
     1690        {
     1691          d = std::sqrt(d2) ;
     1692          s = -b + d;    // second root
     1693
     1694          if ((s >= 0) && (s < snxt))
     1695          {
     1696            xi    = p.x() + s*v.x() ;
     1697            yi    = p.y() + s*v.y() ;
     1698            zi    = p.z() + s*v.z() ;
     1699            rhoi2 = xi*xi + yi*yi   ;
     1700            radi2 = rhoi2 + zi*zi   ;
     1701
     1702            if ( (radi2 <= tolORMax2)
     1703              && (radi2 >= tolORMin2)
     1704              && (zi*(eTheta - halfpi) <= 0) )
     1705            {
     1706              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
     1707              {
     1708                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
     1709                if ( cosPsi >= cosHDPhiOT )
     1710                {
     1711                  snxt=s;
     1712                }
     1713              }
     1714              else
     1715              {
     1716                snxt = s ;
     1717              }
    18261718            }
    18271719          }
    … …  
    18441736{
    18451737  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
    1846   G4double rho2,rad,rho;
    1847   G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi;
     1738  G4double rho2,rds,rho;
     1739  G4double cosPsi;
    18481740  G4double pTheta,dTheta1,dTheta2;
    18491741  rho2=p.x()*p.x()+p.y()*p.y();
    1850   rad=std::sqrt(rho2+p.z()*p.z());
     1742  rds=std::sqrt(rho2+p.z()*p.z());
    18511743  rho=std::sqrt(rho2);
    18521744
    … …  
    18561748  if (fRmin)
    18571749  {
    1858     safeRMin=fRmin-rad;
    1859     safeRMax=rad-fRmax;
     1750    safeRMin=fRmin-rds;
     1751    safeRMax=rds-fRmax;
    18601752    if (safeRMin>safeRMax)
    18611753    {
    … …  
    18691761  else
    18701762  {
    1871     safe=rad-fRmax;
     1763    safe=rds-fRmax;
    18721764  }
    18731765
    … …  
    18751767  // Distance to phi extent
    18761768  //
    1877   if (fDPhi<twopi&&rho)
    1878   {
    1879     phiC=fSPhi+fDPhi*0.5;
    1880     cosPhiC=std::cos(phiC);
    1881     sinPhiC=std::sin(phiC);
    1882 
     1769  if (!fFullPhiSphere && rho)
     1770  {
    18831771    // Psi=angle from central phi to point
    18841772    //
    1885     cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
    1886     if (cosPsi<std::cos(fDPhi*0.5))
     1773    cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
     1774    if (cosPsi<std::cos(hDPhi))
    18871775    {
    18881776      // Point lies outside phi range
    18891777      //
    1890       if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
    1891       {
    1892         safePhi=std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
     1778      if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
     1779      {
     1780        safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
    18931781      }
    18941782      else
    18951783      {
    1896         ePhi=fSPhi+fDPhi;
    1897         safePhi=std::fabs(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
    1898       }
    1899       if (safePhi>safe) safe=safePhi;
     1784        safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
     1785      }
     1786      if (safePhi>safe)  { safe=safePhi; }
    19001787    }
    19011788  }
    … …  
    19031790  // Distance to Theta extent
    19041791  //   
    1905   if ((rad!=0.0) && (fDTheta<pi))
    1906   {
    1907     pTheta=std::acos(p.z()/rad);
    1908     if (pTheta<0) pTheta+=pi;
     1792  if ((rds!=0.0) && (!fFullThetaSphere))
     1793  {
     1794    pTheta=std::acos(p.z()/rds);
     1795    if (pTheta<0)  { pTheta+=pi; }
    19091796    dTheta1=fSTheta-pTheta;
    1910     dTheta2=pTheta-(fSTheta+fDTheta);
     1797    dTheta2=pTheta-eTheta;
    19111798    if (dTheta1>dTheta2)
    19121799    {
    19131800      if (dTheta1>=0)             // WHY ???????????
    19141801      {
    1915         safeTheta=rad*std::sin(dTheta1);
     1802        safeTheta=rds*std::sin(dTheta1);
    19161803        if (safe<=safeTheta)
    19171804        {
    … …  
    19241811      if (dTheta2>=0)
    19251812      {
    1926         safeTheta=rad*std::sin(dTheta2);
     1813        safeTheta=rds*std::sin(dTheta2);
    19271814        if (safe<=safeTheta)
    19281815        {
    … …  
    19331820  }
    19341821
    1935   if (safe<0) safe=0;
     1822  if (safe<0)  { safe=0; }
    19361823  return safe;
    19371824}
    … …  
    19391826/////////////////////////////////////////////////////////////////////
    19401827//
    1941 // Calculate distance to surface of shape from `inside', allowing for tolerance
     1828// Calculate distance to surface of shape from 'inside', allowing for tolerance
    19421829// - Only Calc rmax intersection if no valid rmin intersection
    19431830
    … …  
    19521839  ESide side=kNull,sidephi=kNull,sidetheta=kNull; 
    19531840
     1841  static const G4double halfCarTolerance = kCarTolerance*0.5;
     1842  static const G4double halfAngTolerance = kAngTolerance*0.5;
     1843  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
     1844  const G4double halfRminTolerance = fRminTolerance*0.5;
     1845  const G4double Rmax_plus  = fRmax + halfRmaxTolerance;
     1846  const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
    19541847  G4double t1,t2;
    19551848  G4double b,c,d;
    … …  
    19571850  // Variables for phi intersection:
    19581851
    1959   G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
    1960   G4double cPhi,sinCPhi,cosCPhi;
    19611852  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
    19621853   
    … …  
    19661857  G4double xi,yi,zi;      // Intersection point
    19671858
    1968   // G4double Comp; // Phi intersection
    1969 
    1970   G4bool segPhi;        // Phi flag and precalcs
    1971   G4double hDPhi,hDPhiOT,hDPhiIT;
    1972   G4double cosHDPhiOT,cosHDPhiIT;
    1973 
    1974   G4bool segTheta;                             // Theta flag and precals
    1975   G4double tanSTheta=0.,tanETheta=0., rhoSecTheta;
    1976   G4double tanSTheta2=0.,tanETheta2=0.;
     1859  // Theta precals
     1860  //
     1861  G4double rhoSecTheta;
    19771862  G4double dist2STheta, dist2ETheta, distTheta;
    19781863  G4double d2,s;
    19791864
    19801865  // General Precalcs
    1981 
     1866  //
    19821867  rho2 = p.x()*p.x()+p.y()*p.y();
    19831868  rad2 = rho2+p.z()*p.z();
    1984   //  G4double rad=std::sqrt(rad2);
    19851869
    19861870  pTheta = std::atan2(std::sqrt(rho2),p.z());
    … …  
    19891873  pDotV3d = pDotV2d+p.z()*v.z();
    19901874
    1991   // Set phi divided flag and precalcs
    1992 
    1993   if( fDPhi < twopi )
    1994   {
    1995     segPhi=true;
    1996     hDPhi=0.5*fDPhi;    // half delta phi
    1997     cPhi=fSPhi+hDPhi;;
    1998     hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
    1999     hDPhiIT=hDPhi-0.5*kAngTolerance;
    2000     sinCPhi=std::sin(cPhi);
    2001     cosCPhi=std::cos(cPhi);
    2002     cosHDPhiOT=std::cos(hDPhiOT);
    2003     cosHDPhiIT=std::cos(hDPhiIT);
    2004   }
    2005   else
    2006   {
    2007     segPhi=false;
    2008   }
    2009 
    20101875  // Theta precalcs
    20111876   
    2012   if ( fDTheta < pi )
    2013   {
    2014     segTheta  = true;
    2015     tolSTheta = fSTheta - kAngTolerance*0.5;
    2016     tolETheta = fSTheta + fDTheta + kAngTolerance*0.5;
    2017   }
    2018   else segTheta = false;
    2019 
     1877  if ( !fFullThetaSphere )
     1878  {
     1879    tolSTheta = fSTheta - halfAngTolerance;
     1880    tolETheta = eTheta + halfAngTolerance;
     1881  }
    20201882   
    20211883  // Radial Intersections from G4Sphere::DistanceToIn
    … …  
    20341896  //
    20351897  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
    2036   //
    2037   // const G4double  fractionTolerance = 1.0e-12;
    2038 
    2039   const G4double  flexRadMaxTolerance =            // kRadTolerance;
    2040     std::max(kRadTolerance, fEpsilon * fRmax);
    2041 
    2042   const G4double  Rmax_plus = fRmax + flexRadMaxTolerance*0.5;
    2043 
    2044   const G4double  flexRadMinTolerance = std::max(kRadTolerance,
    2045                      fEpsilon * fRmin);
    2046 
    2047   const G4double  Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;
    2048 
    2049   if(rad2 <= Rmax_plus*Rmax_plus && rad2 >= Rmin_minus*Rmin_minus)
    2050     //  if(rad <= Rmax_plus && rad >= Rmin_minus)
     1898
     1899  if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
    20511900  {
    20521901    c = rad2 - fRmax*fRmax;
    20531902
    2054     if (c < flexRadMaxTolerance*fRmax)
     1903    if (c < fRmaxTolerance*fRmax)
    20551904    {
    20561905      // Within tolerant Outer radius
    … …  
    20651914      d2 = pDotV3d*pDotV3d - c;
    20661915
    2067       if( (c >- flexRadMaxTolerance*fRmax)       // on tolerant surface
    2068        && ((pDotV3d >=0) || (d2 < 0)) )          // leaving outside from Rmax
    2069                                                  // not re-entering
     1916      if( (c >- fRmaxTolerance*fRmax)       // on tolerant surface
     1917       && ((pDotV3d >=0) || (d2 < 0)) )     // leaving outside from Rmax
     1918                                            // not re-entering
    20701919      {
    20711920        if(calcNorm)
    … …  
    20921941      d2 = pDotV3d*pDotV3d - c;
    20931942
    2094       if ( c >- flexRadMinTolerance*fRmin ) // 2.0 * (0.5*kRadTolerance) * fRmin
    2095       {
    2096         if( c < flexRadMinTolerance*fRmin &&
    2097             d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
    2098         {
    2099           if(calcNorm)  *validNorm = false ;   // Rmin surface is concave         
    2100                        return snxt = 0 ;
     1943      if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
     1944      {
     1945        if ( (c < fRminTolerance*fRmin)              // leaving from Rmin
     1946          && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
     1947        {
     1948          if(calcNorm)  { *validNorm = false; }  // Rmin surface is concave         
     1949          return snxt = 0 ;
    21011950        }
    21021951        else
    … …  
    21191968  // Theta segment intersection
    21201969
    2121   if (segTheta)
     1970  if ( !fFullThetaSphere )
    21221971  {
    21231972    // Intersection with theta surfaces
    … …  
    21431992    //
    21441993 
    2145     /* ////////////////////////////////////////////////////////
    2146 
    2147     tanSTheta=std::tan(fSTheta);
    2148     tanSTheta2=tanSTheta*tanSTheta;
    2149     tanETheta=std::tan(fSTheta+fDTheta);
    2150     tanETheta2=tanETheta*tanETheta;
    2151      
    2152     if (fSTheta)
    2153     {
    2154       dist2STheta=rho2-p.z()*p.z()*tanSTheta2;
    2155     }
    2156     else
    2157     {
    2158       dist2STheta = kInfinity;
    2159     }
    2160     if (fSTheta + fDTheta < pi)
    2161     {
    2162       dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
    2163     }
    2164     else
    2165     {
    2166       dist2ETheta = kInfinity ;
    2167     }
    2168     if (pTheta > tolSTheta && pTheta < tolETheta)   // Inside theta 
    2169     {
    2170       // In tolerance of STheta and possible leaving out to small thetas N-
    2171 
    2172       if(pTheta < tolSTheta + kAngTolerance  && fSTheta > kAngTolerance) 
    2173       {
    2174         t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta
    2175 
    2176         if( fSTheta < pi*0.5 && t2 < 0)
    2177         {
    2178           if(calcNorm) *validNorm = false ;
    2179           return snxt = 0 ;
    2180         }
    2181         else if(fSTheta > pi*0.5 && t2 >= 0)
    2182         {
    2183           if(calcNorm)
    2184           {
    2185             rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)) ;
    2186             *validNorm = true ;
    2187             *n = G4ThreeVector(-p.x()/rhoSecTheta,   // N-
    2188                                -p.y()/rhoSecTheta,
    2189                                tanSTheta/std::sqrt(1+tanSTheta2) ) ;
    2190           }
    2191           return snxt = 0 ;
    2192         }
    2193         else if( fSTheta == pi*0.5 && v.z() > 0)
    2194         {
    2195           if(calcNorm)
    2196           {
    2197             *validNorm = true ;
    2198             *n = G4ThreeVector(0,0,1) ;
    2199           }
    2200           return snxt = 0 ;
    2201         }
    2202       }
    2203 
    2204       // In tolerance of ETheta and possible leaving out to larger thetas N+
    2205 
    2206       if ( (pTheta  > tolETheta - kAngTolerance)
    2207         && (( fSTheta + fDTheta) < pi - kAngTolerance) ) 
    2208       {
    2209         t2=pDotV2d-p.z()*v.z()*tanETheta2 ;
    2210         if((fSTheta+fDTheta)>pi*0.5 && t2<0)
    2211         {
    2212           if(calcNorm) *validNorm = false ;
    2213           return snxt = 0 ;
    2214         }
    2215         else if( (fSTheta+fDTheta) < pi*0.5 && t2 >= 0 )
    2216         {
    2217           if(calcNorm)
    2218           {
    2219             rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)) ;
    2220             *validNorm = true ;
    2221             *n = G4ThreeVector( p.x()/rhoSecTheta,  // N+
    2222                                 p.y()/rhoSecTheta,
    2223                                 -tanETheta/std::sqrt(1+tanETheta2)  ) ;
    2224           }
    2225           return snxt = 0 ;
    2226         }
    2227         else if( ( fSTheta+fDTheta) == pi*0.5 && v.z() < 0 )
    2228         {
    2229           if(calcNorm)
    2230           {
    2231             *validNorm = true ;
    2232             *n = G4ThreeVector(0,0,-1) ;
    2233           }
    2234           return snxt = 0 ;
    2235         }
    2236       }
    2237       if( fSTheta > 0 )
    2238       {       
    2239         // First root of fSTheta cone, second if first root -ve
    2240 
    2241         t1 = 1-v.z()*v.z()*(1+tanSTheta2);
    2242         t2 = pDotV2d-p.z()*v.z()*tanSTheta2;
    2243        
    2244         b  = t2/t1;
    2245         c  = dist2STheta/t1;
    2246         d2 = b*b - c ;
    2247 
    2248         if ( d2 >= 0 )
    2249         {
    2250           d = std::sqrt(d2) ;
    2251           s = -b - d ;    // First root
    2252 
    2253           if ( s < 0 )
    2254           {
    2255             s = -b + d ;    // Second root
    2256           }
    2257           if (s > flexRadMaxTolerance*0.5 )   // && s<sr)
    2258           {
    2259             // check against double cone solution
    2260             zi=p.z()+s*v.z();
    2261             if (fSTheta<pi*0.5 && zi<0)
    2262             {
    2263               s = kInfinity ;  // wrong cone
    2264             }
    2265             if (fSTheta>pi*0.5 && zi>0)
    2266             {
    2267               s = kInfinity ;  // wrong cone
    2268             }
    2269             stheta = s ;
    2270             sidetheta = kSTheta ;
    2271           }
    2272         }
    2273       }
    2274 
    2275       // Possible intersection with ETheta cone 
    2276      
    2277       if (fSTheta + fDTheta < pi)
    2278       {
    2279         t1 = 1-v.z()*v.z()*(1+tanETheta2);
    2280         t2 = pDotV2d-p.z()*v.z()*tanETheta2;       
    2281         b  = t2/t1;
    2282         c  = dist2ETheta/t1;
    2283         d2 = b*b-c ;
    2284 
    2285         if ( d2 >= 0 )
    2286         {
    2287           d = std::sqrt(d2);
    2288           s = -b - d ;          // First root
    2289 
    2290           if ( s < 0 )
    2291           {
    2292             s=-b+d;    // Second root
    2293           }
    2294           if (s > flexRadMaxTolerance*0.5 && s < stheta )
    2295           {
    2296             // check against double cone solution
    2297             zi=p.z()+s*v.z();
    2298             if (fSTheta+fDTheta<pi*0.5 && zi<0)
    2299             {
    2300               s = kInfinity ;  // wrong cone
    2301             }
    2302             if (fSTheta+fDTheta>pi*0.5 && zi>0)
    2303             {
    2304               s = kInfinity ;  // wrong cone
    2305             }
    2306           }
    2307           if (s < stheta)
    2308           {
    2309             stheta = s ;
    2310             sidetheta = kETheta ;
    2311           }
    2312         }
    2313       }
    2314     } 
    2315     */  ////////////////////////////////////////////////////////////
    2316 
    23171994    if(fSTheta) // intersection with first cons
    23181995    {
    2319 
    2320       tanSTheta = std::tan(fSTheta);
    2321 
    23221996      if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
    23231997      {
    23241998        if( v.z() > 0. )
    23251999        {
    2326           if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
     2000          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
    23272001          {
    23282002            if(calcNorm)
    … …  
    23332007            return snxt = 0 ;
    23342008          } 
    2335           // s = -p.z()/v.z();
    23362009          stheta    = -p.z()/v.z();
    23372010          sidetheta = kSTheta;
    … …  
    23402013      else // kons is not plane
    23412014      {
    2342         tanSTheta2  = tanSTheta*tanSTheta;
    23432015        t1          = 1-v.z()*v.z()*(1+tanSTheta2);
    23442016        t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
    2345         dist2STheta = rho2-p.z()*p.z()*tanSTheta2;      // t3
    2346 
    2347         // distTheta = std::sqrt(std::fabs(dist2STheta/(1+tanSTheta2)));
     2017        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;     // t3
     2018
    23482019        distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
    23492020
    2350         if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
    2351         {
     2021        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
     2022        {                                      // v parallel to kons
    23522023          if( v.z() > 0. )
    23532024          {
    2354             if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
    2355             {
    2356               if( fSTheta < halfpi && p.z() > 0. )
    2357               {
    2358                 if( calcNorm ) *validNorm = false;
    2359                               return snxt = 0.;
    2360               }
    2361               else if( fSTheta > halfpi && p.z() <= 0)
     2025            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
     2026            {
     2027              if( (fSTheta < halfpi) && (p.z() > 0.) )
     2028              {
     2029                if( calcNorm )  { *validNorm = false; }
     2030                return snxt = 0.;
     2031              }
     2032              else if( (fSTheta > halfpi) && (p.z() <= 0) )
    23622033              {
    23632034                if( calcNorm )
    … …  
    23772048              }
    23782049            }
    2379             // s = -0.5*dist2STheta/t2;
    2380 
    23812050            stheta    = -0.5*dist2STheta/t2;
    23822051            sidetheta = kSTheta;
    23832052          } 
    2384         }
    2385         else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
    2386         {
    2387           if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
    2388           {
    2389             if( fSTheta > halfpi && t2 >= 0. ) // leave
     2053        }      // 2nd order equation, 1st root of fSTheta cone,
     2054        else   // 2nd if 1st root -ve
     2055        {
     2056          if( std::fabs(distTheta) < halfRmaxTolerance )
     2057          {
     2058            if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
    23902059            {
    23912060              if( calcNorm )
    … …  
    23942063                if (rho2)
    23952064                {
    2396                     rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
     2065                  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
    23972066                   
    2398                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
    2399                                         p.y()/rhoSecTheta,
    2400                                         std::sin(fSTheta)  );
     2067                  *n = G4ThreeVector( p.x()/rhoSecTheta,   
     2068                                      p.y()/rhoSecTheta,
     2069                                      std::sin(fSTheta)  );
    24012070                }
    2402                 else *n = G4ThreeVector(0.,0.,1.);
    2403               }                           
     2071                else  { *n = G4ThreeVector(0.,0.,1.); }
     2072              }
    24042073              return snxt = 0.;
    24052074            }
    2406             else if( fSTheta < halfpi  && t2 < 0. && p.z() >=0. ) // leave
    2407             {
    2408                 if( calcNorm )   *validNorm = false;                                                 
    2409                 return snxt = 0.;
     2075            else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
     2076            {
     2077              if( calcNorm )  { *validNorm = false; }
     2078              return snxt = 0.;
    24102079            }                               
    24112080          }
    … …  
    24222091              s = -b - d;         // First root
    24232092
    2424               if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
    2425                               s < 0.  ||
    2426                   ( s > 0. && p.z() + s*v.z() > 0.)                   )
     2093              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
     2094               ||  (s < 0.)  || ( (s > 0.) && (p.z() + s*v.z() > 0.) )     )
    24272095              {
    24282096                s = -b + d ; // 2nd root
    24292097              }
    2430               if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 
     2098              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) ) 
    24312099              {
    24322100                stheta    = s;
    … …  
    24382106              s = -b - d;         // First root
    24392107
    2440               if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
    2441                               s < 0.                                   ||
    2442                   ( s > 0. && p.z() + s*v.z() < 0.)                      )
     2108              if ( ( (std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.) )
     2109                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() < 0.) )   )
    24432110              {
    24442111                s = -b + d ; // 2nd root
    24452112              }
    2446               if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() >= 0.) 
     2113              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() >= 0.) ) 
    24472114              {
    24482115                stheta    = s;
    … …  
    24542121      }
    24552122    }
    2456     if (fSTheta + fDTheta < pi) // intersection with second cons
    2457     {
    2458 
    2459       tanETheta = std::tan(fSTheta+fDTheta);
    2460 
     2123    if (eTheta < pi) // intersection with second cons
     2124    {
    24612125      if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
    24622126      {
    24632127        if( v.z() < 0. )
    24642128        {
    2465           if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
     2129          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
    24662130          {
    24672131            if(calcNorm)
    … …  
    24742138          s = -p.z()/v.z();
    24752139
    2476           if( s < stheta)
     2140          if( s < stheta )
    24772141          {
    24782142            stheta    = s;
    … …  
    24832147      else // kons is not plane
    24842148      {
    2485         tanETheta2  = tanETheta*tanETheta;
    24862149        t1          = 1-v.z()*v.z()*(1+tanETheta2);
    24872150        t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
    2488         dist2ETheta = rho2-p.z()*p.z()*tanETheta2;      // t3
    2489 
    2490         // distTheta = std::sqrt(std::fabs(dist2ETheta/(1+tanETheta2)));
     2151        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;     // t3
     2152
    24912153        distTheta = std::sqrt(rho2)-p.z()*tanETheta;
    24922154
    2493         if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
    2494         {
     2155        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
     2156        {                                      // v parallel to kons
    24952157          if( v.z() < 0. )
    24962158          {
    2497             if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
    2498             {
    2499               if( fSTheta+fDTheta > halfpi && p.z() < 0. )
    2500               {
    2501                 if( calcNorm ) *validNorm = false;
    2502                               return snxt = 0.;
    2503               }
    2504               else if( fSTheta+fDTheta < halfpi && p.z() >= 0)
     2159            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
     2160            {
     2161              if( (eTheta > halfpi) && (p.z() < 0.) )
     2162              {
     2163                if( calcNorm )  { *validNorm = false; }
     2164                return snxt = 0.;
     2165              }
     2166              else if ( (eTheta < halfpi) && (p.z() >= 0) )
    25052167              {
    25062168                if( calcNorm )
    … …  
    25102172                  {
    25112173                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
    2512                    
    25132174                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
    25142175                                        p.y()/rhoSecTheta,
    2515                                         -std::sin(fSTheta+fDTheta)  );
     2176                                        -sinETheta  );
    25162177                  }
    2517                   else *n = G4ThreeVector(0.,0.,-1.);
     2178                  else  { *n = G4ThreeVector(0.,0.,-1.); }
    25182179                }
    25192180                return snxt = 0.;               
    … …  
    25222183            s = -0.5*dist2ETheta/t2;
    25232184
    2524             if( s < stheta)
     2185            if( s < stheta )
    25252186            {
    25262187              stheta    = s;
    … …  
    25282189            }
    25292190          } 
    2530         }
    2531         else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
    2532         {
    2533           if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
    2534           {
    2535             if( fSTheta+fDTheta < halfpi && t2 >= 0. ) // leave
     2191        }      // 2nd order equation, 1st root of fSTheta cone
     2192        else   // 2nd if 1st root -ve
     2193        {
     2194          if ( std::fabs(distTheta) < halfRmaxTolerance )
     2195          {
     2196            if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
    25362197            {
    25372198              if( calcNorm )
    … …  
    25412202                {
    25422203                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
    2543                    
    25442204                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
    25452205                                        p.y()/rhoSecTheta,
    2546                                         -std::sin(fSTheta+fDTheta)  );
     2206                                        -sinETheta  );
    25472207                }
    25482208                else *n = G4ThreeVector(0.,0.,-1.);
    … …  
    25502210              return snxt = 0.;
    25512211            }
    2552             else if( fSTheta+fDTheta > halfpi  && t2 < 0. && p.z() <=0. ) // leave
    2553             {
    2554                 if( calcNorm )   *validNorm = false;                                                 
    2555                 return snxt = 0.;
     2212            else if ( (eTheta > halfpi)
     2213                   && (t2 < 0.) && (p.z() <=0.) ) // leave
     2214            {
     2215              if( calcNorm )  { *validNorm = false; }
     2216              return snxt = 0.;
    25562217            }                               
    25572218          }
    … …  
    25642225            d = std::sqrt(d2);
    25652226
    2566             if( fSTheta+fDTheta < halfpi )
     2227            if( eTheta < halfpi )
    25672228            {
    25682229              s = -b - d;         // First root
    25692230
    2570               if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
    2571                               s < 0. )
     2231              if( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
     2232               || (s < 0.) )
    25722233              {
    25732234                s = -b + d ; // 2nd root
    25742235              }
    2575               if( s >  flexRadMaxTolerance*0.5 ) 
     2236              if( s > halfRmaxTolerance ) 
    25762237              {
    25772238                if( s < stheta )
    … …  
    25862247              s = -b - d;         // First root
    25872248
    2588               if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
    2589                               s < 0.                                   ||
    2590                   ( s > 0. && p.z() + s*v.z() > 0.)                      )
     2249              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.))
     2250                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() > 0.) ) )
    25912251              {
    25922252                s = -b + d ; // 2nd root
    25932253              }
    2594               if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 
     2254              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) ) 
    25952255              {
    25962256                if( s < stheta )
    … …  
    26102270  // Phi Intersection
    26112271   
    2612   if ( fDPhi < twopi)
    2613   {
    2614     sinSPhi=std::sin(fSPhi);
    2615     cosSPhi=std::cos(fSPhi);
    2616     ePhi=fSPhi+fDPhi;
    2617     sinEPhi=std::sin(ePhi);
    2618     cosEPhi=std::cos(ePhi);
    2619     cPhi=fSPhi+fDPhi*0.5;
    2620     sinCPhi=std::sin(cPhi);
    2621     cosCPhi=std::cos(cPhi);
    2622 
    2623     if ( p.x()||p.y() ) // Check if on z axis (rho not needed later)
     2272  if ( !fFullPhiSphere )
     2273  {
     2274    if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
    26242275    {
    26252276      // pDist -ve when inside
    … …  
    26342285      sidephi = kNull ;
    26352286
    2636       if ( pDistS <= 0 && pDistE <= 0 )
     2287      if ( (pDistS <= 0) && (pDistE <= 0) )
    26372288      {
    26382289        // Inside both phi *full* planes
    … …  
    26442295          yi   = p.y()+sphi*v.y() ;
    26452296
    2646           // Check intersecting with correct half-plane
    2647           // (if not -> no intersect)
    2648 
    2649           if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
     2297          // Check intersection with correct half-plane (if not -> no intersect)
     2298          //
     2299          if( (std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance) )
     2300          {
     2301            vphi = std::atan2(v.y(),v.x());
     2302            sidephi = kSPhi;
     2303            if ( ( (fSPhi-halfAngTolerance) <= vphi)
     2304              && ( (ePhi+halfAngTolerance)  >= vphi) )
     2305            {
     2306              sphi = kInfinity;
     2307            }
     2308          }
     2309          else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
    26502310          {
    26512311            sphi=kInfinity;
    … …  
    26542314          {
    26552315            sidephi = kSPhi ;
    2656             if ( pDistS > -0.5*kCarTolerance) sphi =0 ; // Leave by sphi
    2657           }
    2658         }
    2659         else sphi = kInfinity ;
     2316            if ( pDistS > -halfCarTolerance)  { sphi = 0; } // Leave by sphi
     2317          }
     2318        }
     2319        else  { sphi = kInfinity; }
    26602320
    26612321        if ( compE < 0 )
    … …  
    26672327            yi = p.y()+sphi2*v.y() ;
    26682328
    2669             // Check intersecting with correct half-plane
    2670  
    2671             if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
     2329            // Check intersection with correct half-plane
     2330            //
     2331            if ((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
     2332            {
     2333              // Leaving via ending phi
     2334              //
     2335              vphi = std::atan2(v.y(),v.x()) ;
     2336               
     2337              if( !((fSPhi-halfAngTolerance <= vphi)
     2338                  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
     2339              {
     2340                sidephi = kEPhi;
     2341                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
     2342                else                                { sphi = 0.0;   }
     2343              }
     2344            }
     2345            else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
    26722346            {
    26732347              sidephi = kEPhi ;
    2674               if ( pDistE <= -0.5*kCarTolerance )
     2348              if ( pDistE <= -halfCarTolerance )
    26752349              {
    26762350                sphi=sphi2;
    … …  
    26842358        }       
    26852359      }
    2686       else if ( pDistS >= 0 && pDistE >= 0 ) // Outside both *full* phi planes
     2360      else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
    26872361      {
    26882362        if ( pDistS <= pDistE )
    … …  
    26962370        if ( fDPhi > pi )
    26972371        {
    2698           if ( compS < 0 && compE < 0 ) sphi = 0 ;
    2699           else                          sphi = kInfinity ;
     2372          if ( (compS < 0) && (compE < 0) )  { sphi = 0; }
     2373          else                               { sphi = kInfinity; }
    27002374        }
    27012375        else
    … …  
    27042378          // will remain inside
    27052379
    2706           if ( compS >= 0 && compE >= 0 )
    2707           {
    2708             sphi=kInfinity;
    2709           }
    2710           else
    2711           {
    2712             sphi=0;
    2713           }
     2380          if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
     2381          else                                { sphi = 0; }
    27142382        }   
    27152383      }
    2716       else if ( pDistS > 0 && pDistE < 0 )
     2384      else if ( (pDistS > 0) && (pDistE < 0) )
    27172385      {
    27182386        // Outside full starting plane, inside full ending plane
    … …  
    27292397            // (if not -> not leaving phi extent)
    27302398            //
    2731             if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
     2399            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
     2400            {
     2401              vphi = std::atan2(v.y(),v.x());
     2402              sidephi = kSPhi;
     2403              if ( ( (fSPhi-halfAngTolerance) <= vphi)
     2404                && ( (ePhi+halfAngTolerance)  >= vphi) )
     2405              {
     2406                sphi = kInfinity;
     2407              }
     2408            }
     2409            else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
    27322410            {
    27332411              sphi = kInfinity ;
    … …  
    27362414            {
    27372415              sidephi = kEPhi ;
    2738               if ( pDistE > -0.5*kCarTolerance ) sphi = 0. ;
     2416              if ( pDistE > -halfCarTolerance )  { sphi = 0.; }
    27392417            }
    27402418          }
    … …  
    27572435              // (if not -> remain in extent)
    27582436              //
    2759               if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
     2437              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
     2438              {
     2439                vphi = std::atan2(v.y(),v.x());
     2440                sidephi = kSPhi;
     2441                if ( ( (fSPhi-halfAngTolerance) <= vphi)
     2442                  && ( (ePhi+halfAngTolerance)  >= vphi) )
     2443                {
     2444                  sphi = kInfinity;
     2445                }
     2446              }
     2447              else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
    27602448              {
    27612449                sphi=kInfinity;
    … …  
    27872475            xi=p.x()+sphi*v.x();
    27882476            yi=p.y()+sphi*v.y();
    2789 
     2477           
    27902478            // Check intersection in correct half-plane
    27912479            // (if not -> not leaving phi extent)
    27922480            //
    2793             if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
     2481            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
     2482            {
     2483              vphi = std::atan2(v.y(),v.x()) ;
     2484              sidephi = kSPhi;
     2485              if ( ( (fSPhi-halfAngTolerance) <= vphi)
     2486                && ( (ePhi+halfAngTolerance)  >= vphi) )
     2487              {
     2488              sphi = kInfinity;
     2489              }
     2490            }
     2491            else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
    27942492            {
    27952493              sphi = kInfinity ;
    27962494            }
    2797            else  // Leaving via Starting phi
    2798            {
     2495            else  // Leaving via Starting phi
     2496            {
    27992497              sidephi = kSPhi ;
    2800              if ( pDistS > -0.5*kCarTolerance ) sphi = 0 ;
     2498              if ( pDistS > -halfCarTolerance )  { sphi = 0; }
    28012499            }
    28022500          }
    … …  
    28152513              xi   = p.x()+sphi*v.x() ;
    28162514              yi   = p.y()+sphi*v.y() ;
    2817 
     2515             
    28182516              // Check intersection in correct half-plane
    28192517              // (if not -> remain in extent)
    28202518              //
    2821               if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
     2519              if((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
     2520              {
     2521                vphi = std::atan2(v.y(),v.x()) ;
     2522                sidephi = kSPhi;
     2523                if ( ( (fSPhi-halfAngTolerance) <= vphi)
     2524                  && ( (ePhi+halfAngTolerance)  >= vphi) )
     2525                {
     2526                  sphi = kInfinity;
     2527                }
     2528              }
     2529              else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
    28222530              {
    28232531                sphi = kInfinity ;
    … …  
    28492557      {
    28502558        vphi = std::atan2(v.y(),v.x()) ;
    2851         if ( fSPhi < vphi && vphi < fSPhi + fDPhi )
    2852         {
    2853           sphi=kInfinity;
     2559        if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
     2560        {
     2561          sphi = kInfinity;
    28542562        }
    28552563        else
    … …  
    28592567        }
    28602568      }
    2861       else  // travel along z - no phi intersaction
     2569      else  // travel along z - no phi intersection
    28622570      {
    28632571        sphi = kInfinity ;
    … …  
    28952603        if ( fDPhi <= pi )     // Normal to Phi-
    28962604        {
    2897           *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
     2605          *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
    28982606          *validNorm=true;
    28992607        }
    2900         else *validNorm=false;
     2608        else  { *validNorm=false; }
    29012609        break ;
    29022610
    … …  
    29042612        if ( fDPhi <= pi )      // Normal to Phi+
    29052613        {
    2906           *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
     2614          *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
    29072615          *validNorm=true;
    29082616        }
    2909         else *validNorm=false;
     2617        else  { *validNorm=false; }
    29102618        break;
    29112619
    … …  
    29202628          xi = p.x() + snxt*v.x();
    29212629          yi = p.y() + snxt*v.y();
    2922           rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2));
    2923           *n = G4ThreeVector( xi/rhoSecTheta,   // N-
    2924                               yi/rhoSecTheta,
    2925                               -tanSTheta/std::sqrt(1+tanSTheta2));
     2630          rho2=xi*xi+yi*yi;
     2631          if (rho2)
     2632          {
     2633            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
     2634            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
     2635                               -tanSTheta/std::sqrt(1+tanSTheta2));
     2636          }
     2637          else
     2638          {
     2639            *n = G4ThreeVector(0.,0.,1.);
     2640          }
    29262641          *validNorm=true;
    29272642        }
    2928         else *validNorm=false;  // Concave STheta cone
     2643        else  { *validNorm=false; }  // Concave STheta cone
    29292644        break;
    29302645
    29312646      case kETheta:
    2932         if( ( fSTheta + fDTheta ) == halfpi )
     2647        if( eTheta == halfpi )
    29332648        {
    29342649          *n         = G4ThreeVector(0.,0.,-1.);
    29352650          *validNorm = true;
    29362651        }
    2937         else if ( ( fSTheta + fDTheta ) < halfpi)
     2652        else if ( eTheta < halfpi )
    29382653        {
    29392654          xi=p.x()+snxt*v.x();
    29402655          yi=p.y()+snxt*v.y();
    2941           rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2));
    2942           *n = G4ThreeVector( xi/rhoSecTheta,   // N+
    2943                               yi/rhoSecTheta,
    2944                               -tanETheta/std::sqrt(1+tanETheta2) );
     2656          rho2=xi*xi+yi*yi;
     2657          if (rho2)
     2658          {
     2659            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
     2660            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
     2661                               -tanETheta/std::sqrt(1+tanETheta2) );
     2662          }
     2663          else
     2664          {
     2665            *n = G4ThreeVector(0.,0.,-1.);
     2666          }
    29452667          *validNorm=true;
    29462668        }
    2947         else *validNorm=false;   // Concave ETheta cone
     2669        else  { *validNorm=false; }   // Concave ETheta cone
    29482670        break;
    29492671
    … …  
    29952717/////////////////////////////////////////////////////////////////////////
    29962718//
    2997 // Calcluate distance (<=actual) to closest surface of shape from inside
     2719// Calculate distance (<=actual) to closest surface of shape from inside
    29982720
    29992721G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
    30002722{
    30012723  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
    3002   G4double rho2,rad,rho;
    3003   G4double phiC,cosPhiC,sinPhiC,ePhi;
     2724  G4double rho2,rds,rho;
    30042725  G4double pTheta,dTheta1,dTheta2;
    30052726  rho2=p.x()*p.x()+p.y()*p.y();
    3006   rad=std::sqrt(rho2+p.z()*p.z());
     2727  rds=std::sqrt(rho2+p.z()*p.z());
    30072728  rho=std::sqrt(rho2);
    30082729
    … …  
    30272748  if (fRmin)
    30282749  {
    3029     safeRMin=rad-fRmin;
    3030     safeRMax=fRmax-rad;
     2750    safeRMin=rds-fRmin;
     2751    safeRMax=fRmax-rds;
    30312752    if (safeRMin<safeRMax)
    30322753    {
    … …  
    30402761  else
    30412762  {
    3042     safe=fRmax-rad;
     2763    safe=fRmax-rds;
    30432764  }
    30442765
    … …  
    30462767  // Distance to phi extent
    30472768  //
    3048   if (fDPhi<twopi && rho)
    3049   {
    3050     phiC=fSPhi+fDPhi*0.5;
    3051     cosPhiC=std::cos(phiC);
    3052     sinPhiC=std::sin(phiC);
    3053     if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
    3054     {
    3055       safePhi=-(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
     2769  if (!fFullPhiSphere && rho)
     2770  {
     2771    if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
     2772    {
     2773      safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
    30562774    }
    30572775    else
    30582776    {
    3059       ePhi=fSPhi+fDPhi;
    3060       safePhi=(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
    3061     }
    3062     if (safePhi<safe) safe=safePhi;
     2777      safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
     2778    }
     2779    if (safePhi<safe)  { safe=safePhi; }
    30632780  }
    30642781
    … …  
    30662783  // Distance to Theta extent
    30672784  //   
    3068   if (rad)
    3069   {
    3070     pTheta=std::acos(p.z()/rad);
    3071     if (pTheta<0) pTheta+=pi;
     2785  if (rds)
     2786  {
     2787    pTheta=std::acos(p.z()/rds);
     2788    if (pTheta<0)  { pTheta+=pi; }
    30722789    dTheta1=pTheta-fSTheta;
    3073     dTheta2=(fSTheta+fDTheta)-pTheta;
    3074     if (dTheta1<dTheta2)
    3075     {
    3076       safeTheta=rad*std::sin(dTheta1);
    3077       if (safe>safeTheta)
    3078       {
    3079         safe=safeTheta;
    3080       }
    3081     }
    3082     else
    3083     {
    3084       safeTheta=rad*std::sin(dTheta2);
    3085       if (safe>safeTheta)
    3086       {
    3087         safe=safeTheta;
    3088       }
    3089     }
    3090   }
    3091 
    3092   if (safe<0) safe=0;
    3093     return safe;
     2790    dTheta2=eTheta-pTheta;
     2791    if (dTheta1<dTheta2)  { safeTheta=rds*std::sin(dTheta1); }
     2792    else                  { safeTheta=rds*std::sin(dTheta2); }
     2793    if (safe>safeTheta)   { safe=safeTheta; }
     2794  }
     2795
     2796  if (safe<0)  { safe=0; }
     2797  return safe;
    30942798}
    30952799
    … …  
    31192823  // Phi cross sections
    31202824   
    3121   noPhiCrossSections=G4int (fDPhi/kMeshAngleDefault)+1;
     2825  noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
    31222826   
    31232827  if (noPhiCrossSections<kMinMeshSections)
    … …  
    31342838  // on the x axis. Will give better extent calculations when not rotated.
    31352839   
    3136   if (fDPhi==pi*2.0 && fSPhi==0)
     2840  if (fFullPhiSphere)
    31372841  {
    31382842    sAnglePhi = -meshAnglePhi*0.5;
    … …  
    31602864  // on the z axis. Will give better extent calculations when not rotated.
    31612865   
    3162   if (fDTheta==pi && fSTheta==0)
     2866  if (fFullThetaSphere)
    31632867  {
    31642868    startTheta = -meshTheta*0.5;
    … …  
    32232927  }
    32242928
    3225   delete[] cosCrossTheta;
    3226   delete[] sinCrossTheta;
     2929  delete [] cosCrossTheta;
     2930  delete [] sinCrossTheta;
    32272931
    32282932  return vertices;
    … …  
    32692973  G4double height1, height2, slant1, slant2, costheta, sintheta,theta,rRand;
    32702974
    3271   height1 = (fRmax-fRmin)*std::cos(fSTheta);
    3272   height2 = (fRmax-fRmin)*std::cos(fSTheta+fDTheta);
    3273   slant1  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta))
    3274                       + height1*height1);
    3275   slant2  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta+fDTheta))
    3276                       + height2*height2);
     2975  height1 = (fRmax-fRmin)*cosSTheta;
     2976  height2 = (fRmax-fRmin)*cosETheta;
     2977  slant1  = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
     2978  slant2  = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
    32772979  rRand   = RandFlat::shoot(fRmin,fRmax);
    32782980 
    3279   aOne = fRmax*fRmax*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
    3280   aTwo = fRmin*fRmin*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
    3281   aThr = fDPhi*((fRmax + fRmin)*std::sin(fSTheta))*slant1;
    3282   aFou = fDPhi*((fRmax + fRmin)*std::sin(fSTheta+fDTheta))*slant2;
     2981  aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
     2982  aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
     2983  aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
     2984  aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
    32832985  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
    32842986 
    3285   phi = RandFlat::shoot(fSPhi, fSPhi + fDPhi);
     2987  phi = RandFlat::shoot(fSPhi, ePhi);
    32862988  cosphi = std::cos(phi);
    32872989  sinphi = std::sin(phi);
    3288   theta = RandFlat::shoot(fSTheta,fSTheta+fDTheta);
     2990  theta = RandFlat::shoot(fSTheta,eTheta);
    32892991  costheta = std::cos(theta);
    32902992  sintheta = std::sqrt(1.-sqr(costheta));
    32912993
    3292   if( ((fSPhi==0) && (fDPhi==2.*pi)) || (fDPhi==2.*pi) ) {aFiv = 0;}
    3293   if(fSTheta == 0)  {aThr=0;}
    3294   if(fDTheta + fSTheta == pi) {aFou = 0;}
    3295   if(fSTheta == 0.5*pi) {aThr = pi*(fRmax*fRmax-fRmin*fRmin);}
    3296   if(fSTheta + fDTheta == 0.5*pi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin);}
     2994  if(fFullPhiSphere) { aFiv = 0; }
     2995  if(fSTheta == 0)   { aThr=0; }
     2996  if(eTheta == pi) { aFou = 0; }
     2997  if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
     2998  if(eTheta == halfpi)  { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
    32972999
    32983000  chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
    … …  
    33093011  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
    33103012  {
    3311     if (fSTheta != 0.5*pi)
    3312     {
    3313       zRand = RandFlat::shoot(fRmin*std::cos(fSTheta),fRmax*std::cos(fSTheta));
    3314       return G4ThreeVector(std::tan(fSTheta)*zRand*cosphi,
    3315                            std::tan(fSTheta)*zRand*sinphi,zRand);
     3013    if (fSTheta != halfpi)
     3014    {
     3015      zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
     3016      return G4ThreeVector(tanSTheta*zRand*cosphi,
     3017                           tanSTheta*zRand*sinphi,zRand);
    33163018    }
    33173019    else
    … …  
    33223024  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
    33233025  {
    3324     if(fSTheta + fDTheta != 0.5*pi)
    3325     {
    3326       zRand = RandFlat::shoot(fRmin*std::cos(fSTheta+fDTheta),
    3327                               fRmax*std::cos(fSTheta+fDTheta));
    3328       return G4ThreeVector  (std::tan(fSTheta+fDTheta)*zRand*cosphi,
    3329                              std::tan(fSTheta+fDTheta)*zRand*sinphi,zRand);
     3026    if(eTheta != halfpi)
     3027    {
     3028      zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
     3029      return G4ThreeVector  (tanETheta*zRand*cosphi,
     3030                             tanETheta*zRand*sinphi,zRand);
    33303031    }
    33313032    else
    … …  
    33363037  else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
    33373038  {
    3338     return G4ThreeVector(rRand*sintheta*std::cos(fSPhi),
    3339                          rRand*sintheta*std::sin(fSPhi),rRand*costheta);
     3039    return G4ThreeVector(rRand*sintheta*cosSPhi,
     3040                         rRand*sintheta*sinSPhi,rRand*costheta);
    33403041  }
    33413042  else
    33423043  {
    3343     return G4ThreeVector(rRand*sintheta*std::cos(fSPhi+fDPhi),
    3344                          rRand*sintheta*std::sin(fSPhi+fDPhi),rRand*costheta);
    3345   }
     3044    return G4ThreeVector(rRand*sintheta*cosEPhi,
     3045                         rRand*sintheta*sinEPhi,rRand*costheta);
     3046  }
     3047}
     3048
     3049////////////////////////////////////////////////////////////////////////////////
     3050//
     3051// GetSurfaceArea
     3052
     3053G4double G4Sphere::GetSurfaceArea()
     3054{
     3055  if(fSurfaceArea != 0.) {;}
     3056  else
     3057  {   
     3058    G4double Rsq=fRmax*fRmax;
     3059    G4double rsq=fRmin*fRmin;
     3060         
     3061    fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
     3062    if(!fFullPhiSphere)
     3063    {
     3064      fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
     3065    }
     3066    if(fSTheta >0)
     3067    {
     3068      G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
     3069                              + std::pow(cosSTheta,2));
     3070      if(fDPhi>pi)
     3071      {
     3072        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
     3073      }
     3074      else
     3075      {
     3076        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
     3077      }
     3078    }
     3079    if(eTheta < pi)
     3080    {
     3081      G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
     3082                              + std::pow(cosETheta,2));
     3083      if(fDPhi>pi)
     3084      {
     3085        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
     3086      }
     3087      else
     3088      {
     3089        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
     3090      }
     3091    }
     3092  }
     3093  return fSurfaceArea;
    33463094}
    33473095

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Torus.cc,v 1.63 2007/10/02 09:34:17 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Torus.cc,v 1.65 2009/11/26 10:31:06 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    284284
    285285    G4double theta = std::atan2(ptmp.y(),ptmp.x());
    286 
    287     if (theta < 0)  { theta += twopi; }
    288286   
     287    if ( fSPhi >= 0 )
     288    {
     289      if ( theta < - kAngTolerance*0.5 )  { theta += twopi; }
     290      if ( (std::abs(theta) < kAngTolerance*0.5)
     291        && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
     292      {
     293        theta += twopi ; // 0 <= theta < 2pi
     294      }
     295    }
     296    if ((fSPhi <= -pi )&&(theta>kAngTolerance*0.5)) { theta = theta-twopi; }
     297       
    289298    // We have to verify if this root is inside the region between
    290299    // fSPhi and fSPhi + fDPhi

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trap.cc,v 1.45 2008/04/23 09:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030// class G4Trap

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

    r1058 r1228  
    2626//
    2727// $Id: G4Trd.cc,v 1.34 2006/10/19 15:33:38 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//

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

    r1058 r1228  
    2525//
    2626//
    27 // $Id: G4Tubs.cc,v 1.74 2008/11/06 15:26:53 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     27// $Id: G4Tubs.cc,v 1.79 2009/06/30 10:10:11 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03 $
    2929//
    3030//
    … …  
    9090                      G4double pDz,
    9191                      G4double pSPhi, G4double pDPhi )
    92   : G4CSGSolid(pName)
     92  : G4CSGSolid(pName), fSPhi(0), fDPhi(0)
    9393{
    9494
    … …  
    102102  else
    103103  {
    104     G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
    105            << "        Negative Z half-length ! - "
    106            << pDz << G4endl;
     104    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
     105           << "        Negative Z half-length (" << pDz << ") in solid: "
     106           << GetName() << G4endl;
    107107    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
    108108                "Invalid Z half-length");
    … …  
    115115  else
    116116  {
    117     G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
    118            << "        Invalid values for radii !" << G4endl
     117    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
     118           << "        Invalid values for radii in solid " << GetName()
     119           << G4endl
    119120           << "        pRMin = " << pRMin << ", pRMax = " << pRMax << G4endl;
    120121    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
    … …  
    122123  }
    123124
    124   fPhiFullTube = true;
    125   if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
    126   {
    127     fDPhi=twopi;
    128     fSPhi=0;
    129   }
    130   else
    131   {
    132     fPhiFullTube = false;
    133     if ( pDPhi > 0 )
    134     {
    135       fDPhi = pDPhi;
    136     }
    137     else
    138     {
    139       G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
    140              << "        Negative delta-Phi ! - "
    141              << pDPhi << G4endl;
    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();
     125  // Check angles
     126
     127  CheckPhiAngles(pSPhi, pDPhi);
    162128}
    163129
    … …  
    203169{
    204170
    205   if ( !pTransform.IsRotated() && fDPhi == twopi && fRMin == 0 )
     171  if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
    206172  {
    207173    // Special case handling for unrotated solid tubes
    … …  
    210176    // and compute exact x and y limit for x/y case
    211177     
    212     G4double xoffset, xMin, xMax ;
    213     G4double yoffset, yMin, yMax ;
    214     G4double zoffset, zMin, zMax ;
    215 
    216     G4double diff1, diff2, maxDiff, newMin, newMax ;
    217     G4double xoff1, xoff2, yoff1, yoff2, delta ;
    218 
    219     xoffset = pTransform.NetTranslation().x() ;
    220     xMin = xoffset - fRMax ;
    221     xMax = xoffset + fRMax ;
     178    G4double xoffset, xMin, xMax;
     179    G4double yoffset, yMin, yMax;
     180    G4double zoffset, zMin, zMax;
     181
     182    G4double diff1, diff2, maxDiff, newMin, newMax;
     183    G4double xoff1, xoff2, yoff1, yoff2, delta;
     184
     185    xoffset = pTransform.NetTranslation().x();
     186    xMin = xoffset - fRMax;
     187    xMax = xoffset + fRMax;
    222188
    223189    if (pVoxelLimit.IsXLimited())
    … …  
    230196      else
    231197      {
    232         if ( xMin < pVoxelLimit.GetMinXExtent() )
    233         {
    234           xMin = pVoxelLimit.GetMinXExtent() ;
    235         }
    236         if (xMax > pVoxelLimit.GetMaxXExtent() )
    237         {
    238           xMax = pVoxelLimit.GetMaxXExtent() ;
    239         }
    240       }
    241     }
    242     yoffset = pTransform.NetTranslation().y() ;
    243     yMin    = yoffset - fRMax ;
    244     yMax    = yoffset + fRMax ;
     198        if (xMin < pVoxelLimit.GetMinXExtent())
     199        {
     200          xMin = pVoxelLimit.GetMinXExtent();
     201        }
     202        if (xMax > pVoxelLimit.GetMaxXExtent())
     203        {
     204          xMax = pVoxelLimit.GetMaxXExtent();
     205        }
     206      }
     207    }
     208    yoffset = pTransform.NetTranslation().y();
     209    yMin    = yoffset - fRMax;
     210    yMax    = yoffset + fRMax;
    245211
    246212    if ( pVoxelLimit.IsYLimited() )
    … …  
    249215        || (yMax < pVoxelLimit.GetMinYExtent()) )
    250216      {
    251         return false ;
     217        return false;
    252218      }
    253219      else
    254220      {
    255         if ( yMin < pVoxelLimit.GetMinYExtent() )
    256         {
    257           yMin = pVoxelLimit.GetMinYExtent() ;
    258         }
    259         if ( yMax > pVoxelLimit.GetMaxYExtent() )
     221        if (yMin < pVoxelLimit.GetMinYExtent())
     222        {
     223          yMin = pVoxelLimit.GetMinYExtent();
     224        }
     225        if (yMax > pVoxelLimit.GetMaxYExtent())
    260226        {
    261227          yMax=pVoxelLimit.GetMaxYExtent();
    … …  
    263229      }
    264230    }
    265     zoffset = pTransform.NetTranslation().z() ;
    266     zMin    = zoffset - fDz ;
    267     zMax    = zoffset + fDz ;
     231    zoffset = pTransform.NetTranslation().z();
     232    zMin    = zoffset - fDz;
     233    zMax    = zoffset + fDz;
    268234
    269235    if ( pVoxelLimit.IsZLimited() )
    … …  
    272238        || (zMax < pVoxelLimit.GetMinZExtent()) )
    273239      {
    274         return false ;
     240        return false;
    275241      }
    276242      else
    277243      {
    278         if ( zMin < pVoxelLimit.GetMinZExtent() )
    279         {
    280           zMin = pVoxelLimit.GetMinZExtent() ;
    281         }
    282         if ( zMax > pVoxelLimit.GetMaxZExtent() )
     244        if (zMin < pVoxelLimit.GetMinZExtent())
     245        {
     246          zMin = pVoxelLimit.GetMinZExtent();
     247        }
     248        if (zMax > pVoxelLimit.GetMaxZExtent())
    283249        {
    284250          zMax = pVoxelLimit.GetMaxZExtent();
    … …  
    290256      case kXAxis :
    291257      {
    292         yoff1 = yoffset - yMin ;
    293         yoff2 = yMax    - yoffset ;
     258        yoff1 = yoffset - yMin;
     259        yoff2 = yMax    - yoffset;
    294260
    295261        if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
    296262        {                                   // => no change
    297           pMin = xMin ;
    298           pMax = xMax ;
     263          pMin = xMin;
     264          pMax = xMax;
    299265        }
    300266        else
    … …  
    317283      case kYAxis :
    318284      {
    319         xoff1 = xoffset - xMin ;
    320         xoff2 = xMax - xoffset ;
     285        xoff1 = xoffset - xMin;
     286        xoff2 = xMax - xoffset;
    321287
    322288        if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
    323289        {                                   // => no change
    324           pMin = yMin ;
    325           pMax = yMax ;
     290          pMin = yMin;
     291          pMax = yMax;
    326292        }
    327293        else
    … …  
    334300          delta   = fRMax*fRMax - xoff2*xoff2;
    335301          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
    336           maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
    337           newMin  = yoffset - maxDiff ;
    338           newMax  = yoffset + maxDiff ;
    339           pMin    = (newMin < yMin) ? yMin : newMin ;
    340           pMax     =(newMax > yMax) ? yMax : newMax ;
    341         }
    342         break ;
     302          maxDiff = (diff1 > diff2) ? diff1 : diff2;
     303          newMin  = yoffset - maxDiff;
     304          newMax  = yoffset + maxDiff;
     305          pMin    = (newMin < yMin) ? yMin : newMin;
     306          pMax    = (newMax > yMax) ? yMax : newMax;
     307        }
     308        break;
    343309      }
    344310      case kZAxis:
    345311      {
    346         pMin = zMin ;
    347         pMax = zMax ;
    348         break ;
     312        pMin = zMin;
     313        pMax = zMax;
     314        break;
    349315      }
    350316      default:
    351317        break;
    352318    }
    353     pMin -= kCarTolerance ;
    354     pMax += kCarTolerance ;
    355     return true;   
     319    pMin -= kCarTolerance;
     320    pMax += kCarTolerance;
     321    return true;
    356322  }
    357323  else // Calculate rotated vertex coordinates
    358324  {
    359     G4int i, noEntries, noBetweenSections4 ;
    360     G4bool existsAfterClip = false ;
    361     G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
     325    G4int i, noEntries, noBetweenSections4;
     326    G4bool existsAfterClip = false;
     327    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
     328
     329    pMin =  kInfinity;
     330    pMax = -kInfinity;
     331
     332    noEntries = vertices->size();
     333    noBetweenSections4 = noEntries - 4;
    362334   
    363     pMin = +kInfinity ;
    364     pMax = -kInfinity ;
    365 
    366     noEntries = vertices->size() ;
    367     noBetweenSections4 = noEntries - 4 ;
    368    
    369     for (i = 0 ; i < noEntries ; i += 4 )
    370     {
    371       ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    372     }
    373     for (i = 0 ; i < noBetweenSections4 ; i += 4 )
    374     {
    375       ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    376     }
    377     if ((pMin != kInfinity) || (pMax != -kInfinity) )
    378     {
    379       existsAfterClip = true ;
    380       pMin -= kCarTolerance ; // Add 2*tolerance to avoid precision troubles
    381       pMax += kCarTolerance ;
     335    for ( i = 0 ; i < noEntries ; i += 4 )
     336    {
     337      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
     338    }
     339    for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
     340    {
     341      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
     342    }
     343    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
     344    {
     345      existsAfterClip = true;
     346      pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
     347      pMax += kCarTolerance;
    382348    }
    383349    else
    … …  
    391357             (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
    392358             (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
    393              (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
     359             (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
    394360       
    395361      if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
    396362      {
    397         existsAfterClip = true ;
    398         pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
    399         pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
     363        existsAfterClip = true;
     364        pMin            = pVoxelLimit.GetMinExtent(pAxis);
     365        pMax            = pVoxelLimit.GetMaxExtent(pAxis);
    400366      }
    401367    }
    … …  
    808774  G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
    809775  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
     776  const G4double dRmax = 100.*fRMax;
    810777
    811778  static const G4double halfCarTolerance = 0.5*kCarTolerance;
    … …  
    908875        if (s >= 0)  // If 'forwards'
    909876        {
     877          if ( s>dRmax ) // Avoid rounding errors due to precision issues on
     878          {              // 64 bits systems. Split long distances and recompute
     879            G4double fTerm = s-std::fmod(s,dRmax);
     880            s = fTerm + DistanceToIn(p+fTerm*v,v);
     881          }
    910882          // Check z intersection
    911883          //
    … …  
    1022994          //
    1023995          if(s < 0.0)  { s = 0.0; }
     996          if ( s>dRmax ) // Avoid rounding errors due to precision issues seen
     997          {              // 64 bits systems. Split long distances and recompute
     998            G4double fTerm = s-std::fmod(s,dRmax);
     999            s = fTerm + DistanceToIn(p+fTerm*v,v);
     1000          }
    10241001          zi = p.z() + s*v.z() ;
    10251002          if (std::fabs(zi) <= tolODz)
Note: See TracChangeset for help on using the changeset viewer.