Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Cons.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4Cons.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Cons.cc (Version 4.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 //
                                                   >>  24 // $Id: G4Cons.cc,v 1.20 2002/01/10 15:42:25 gcosmo Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-04-01 $
                                                   >>  26 //
                                                   >>  27 // class G4Cons
                                                   >>  28 //
 26 // Implementation for G4Cons class                 29 // Implementation for G4Cons class
 27 //                                                 30 //
 28 // ~1994    P.Kent: Created, as main part of t <<  31 // History:
 29 // 13.09.96 V.Grichine: Review and final modif <<  32 //
 30 // 03.05.05 V.Grichine: SurfaceNormal(p) accor <<  33 // 05.10.00 V. Grichine bugs fixed in   Distance ToIn(p,v)
 31 // 12.10.09 T.Nikitina: Added to DistanceToIn( <<  34 // 17.08.00 V.Grichine, if one and only one Rmin=0, it'll be 1e3*kRadTolerance 
 32 //                      case of point on surfa <<  35 // 08.08.00 V.Grichine, more stable roots of 2-equation in DistanceToOut(p,v,...) 
 33 // 03.10.16 E.Tcherniaev: use G4BoundingEnvelo <<  36 // 06.03.00 V.Grichine, modifications in DistanceToOut(p,v,...) 
 34 //                      removed CreateRotatedV <<  37 // 18.11.99 V.Grichine side = kNull initialisation in DistanceToOut(p,v,...)
 35 // ------------------------------------------- <<  38 // 28.04.99 V. Grichine bugs fixed in  Distance ToOut(p,v,...) and  
                                                   >>  39 //          Distance ToIn(p,v)
                                                   >>  40 // 09.10.98 V. Grichine modifications in Distance ToOut(p,v,...)
                                                   >>  41 // 13.09.96 V. Grichine: final modifications to commit
                                                   >>  42 // ~1994 P. Kent: main part of geometry functions
 36                                                    43 
 37 #include "G4Cons.hh"                               44 #include "G4Cons.hh"
 38                                                    45 
 39 #if !defined(G4GEOM_USE_UCONS)                 << 
 40                                                << 
 41 #include "G4GeomTools.hh"                      << 
 42 #include "G4VoxelLimits.hh"                        46 #include "G4VoxelLimits.hh"
 43 #include "G4AffineTransform.hh"                    47 #include "G4AffineTransform.hh"
 44 #include "G4BoundingEnvelope.hh"               << 
 45 #include "G4GeometryTolerance.hh"              << 
 46                                                    48 
 47 #include "G4VPVParameterisation.hh"                49 #include "G4VPVParameterisation.hh"
 48                                                    50 
 49 #include "meshdefs.hh"                             51 #include "meshdefs.hh"
 50                                                    52 
 51 #include "Randomize.hh"                        << 
 52                                                << 
 53 #include "G4VGraphicsScene.hh"                     53 #include "G4VGraphicsScene.hh"
                                                   >>  54 #include "G4Polyhedron.hh"
                                                   >>  55 #include "G4NURBS.hh"
                                                   >>  56 #include "G4NURBSbox.hh"
 54                                                    57 
 55 using namespace CLHEP;                         << 
 56                                                << 
 57 //////////////////////////////////////////////     58 ////////////////////////////////////////////////////////////////////////
 58 //                                                 59 //
 59 // Private enums: Not for external use         <<  60 // Private enum: Not for external use - used by distanceToOut
 60                                                << 
 61 namespace                                      << 
 62 {                                              << 
 63   // used by DistanceToOut()                   << 
 64   //                                           << 
 65   enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kP << 
 66                                                << 
 67   // used by ApproxSurfaceNormal()             << 
 68   //                                           << 
 69   enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ} << 
 70 }                                              << 
 71                                                    61 
 72 ////////////////////////////////////////////// <<  62 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
 73 //                                             << 
 74 // constructor - check parameters, convert ang << 
 75 //             - note if pDPhi>2PI then reset  << 
 76                                                    63 
 77 G4Cons::G4Cons( const G4String& pName,         <<  64 // used by normal
 78                       G4double  pRmin1, G4doub << 
 79                       G4double  pRmin2, G4doub << 
 80                       G4double pDz,            << 
 81                       G4double pSPhi, G4double << 
 82   : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2( << 
 83     fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz),  << 
 84 {                                              << 
 85   kRadTolerance = G4GeometryTolerance::GetInst << 
 86   kAngTolerance = G4GeometryTolerance::GetInst << 
 87                                                    65 
 88   halfCarTolerance=kCarTolerance*0.5;          <<  66 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
 89   halfRadTolerance=kRadTolerance*0.5;          << 
 90   halfAngTolerance=kAngTolerance*0.5;          << 
 91                                                << 
 92   // Check z-len                               << 
 93   //                                           << 
 94   if ( pDz < 0 )                               << 
 95   {                                            << 
 96     std::ostringstream message;                << 
 97     message << "Invalid Z half-length for Soli << 
 98             << "        hZ = " << pDz;         << 
 99     G4Exception("G4Cons::G4Cons()", "GeomSolid << 
100                 FatalException, message);      << 
101   }                                            << 
102                                                << 
103   // Check radii                               << 
104   //                                           << 
105   if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || << 
106   {                                            << 
107     std::ostringstream message;                << 
108     message << "Invalid values of radii for So << 
109             << "        pRmin1 = " << pRmin1 < << 
110             << ", pRmax1 = " << pRmax1 << ", p << 
111     G4Exception("G4Cons::G4Cons()", "GeomSolid << 
112                 FatalException, message) ;     << 
113   }                                            << 
114   if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fR << 
115   if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fR << 
116                                                << 
117   // Check angles                              << 
118   //                                           << 
119   CheckPhiAngles(pSPhi, pDPhi);                << 
120 }                                              << 
121                                                    67 
122 //////////////////////////////////////////////     68 ///////////////////////////////////////////////////////////////////////
123 //                                                 69 //
124 // Fake default constructor - sets only member <<  70 // Destructor
125 //                            for usage restri <<  71 
126 //                                             <<  72 G4Cons::~G4Cons()
127 G4Cons::G4Cons( __void__& a )                  << 
128   : G4CSGSolid(a)                              << 
129 {                                                  73 {
                                                   >>  74    ;
130 }                                                  75 }
131                                                    76 
132 ////////////////////////////////////////////// <<  77 //////////////////////////////////////////////////////////////////////////
133 //                                                 78 //
134 // Destructor                                  <<  79 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
                                                   >>  80 //               - note if pDPhi>2PI then reset to 2PI
135                                                    81 
136 G4Cons::~G4Cons() = default;                   <<  82 G4Cons::G4Cons( const G4String& pName,
                                                   >>  83           G4double pRmin1, G4double pRmax1,
                                                   >>  84           G4double pRmin2, G4double pRmax2,
                                                   >>  85           G4double pDz,
                                                   >>  86           G4double pSPhi, G4double pDPhi) : G4CSGSolid(pName)
                                                   >>  87 {
                                                   >>  88 // Check z-len
137                                                    89 
138 ////////////////////////////////////////////// <<  90   if ( pDz > 0 ) fDz = pDz ;
139 //                                             <<  91   else
140 // Copy constructor                            <<  92   {
                                                   >>  93     G4Exception("Error in G4Cons::G4Cons - invalid z half-length");
                                                   >>  94   }
141                                                    95 
142 G4Cons::G4Cons(const G4Cons&) = default;       <<  96 // Check radii
143                                                    97 
144 ////////////////////////////////////////////// <<  98   if ( pRmin1 < pRmax1 && pRmin2 < pRmax2 && pRmin1 >= 0 && pRmin2 >= 0 )
145 //                                             <<  99   {
146 // Assignment operator                         << 
147                                                   100 
148 G4Cons& G4Cons::operator = (const G4Cons& rhs) << 101     fRmin1 = pRmin1 ; 
149 {                                              << 102     fRmax1 = pRmax1 ;
150    // Check assignment to self                 << 103     fRmin2 = pRmin2 ; 
151    //                                          << 104     fRmax2 = pRmax2 ;
152    if (this == &rhs)  { return *this; }        << 105     if( (pRmin1 == 0.0 && pRmin2 > 0.0) ) fRmin1 = 1e3*kRadTolerance ; 
153                                                << 106     if( (pRmin2 == 0.0 && pRmin1 > 0.0) ) fRmin2 = 1e3*kRadTolerance ; 
154    // Copy base class data                     << 107   }
155    //                                          << 108   else  G4Exception("Error in G4Cons::G4Cons - invalid radii") ;
156    G4CSGSolid::operator=(rhs);                 << 
157                                                << 
158    // Copy data                                << 
159    //                                          << 
160    kRadTolerance = rhs.kRadTolerance;          << 
161    kAngTolerance = rhs.kAngTolerance;          << 
162    fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;   << 
163    fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;   << 
164    fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = r << 
165    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh << 
166    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r << 
167    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh << 
168    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh << 
169    fPhiFullCone = rhs.fPhiFullCone;            << 
170    halfCarTolerance = rhs.halfCarTolerance;    << 
171    halfRadTolerance = rhs.halfRadTolerance;    << 
172    halfAngTolerance = rhs.halfAngTolerance;    << 
173                                                   109 
174    return *this;                               << 110 // Check angles
                                                   >> 111 
                                                   >> 112   if ( pDPhi >= 2.0*M_PI )
                                                   >> 113   {
                                                   >> 114        fDPhi=2*M_PI;
                                                   >> 115        fSPhi=0;
                                                   >> 116   }
                                                   >> 117   else
                                                   >> 118   {
                                                   >> 119     if ( pDPhi > 0 ) fDPhi = pDPhi ;
                                                   >> 120     else             G4Exception("Error in G4Cons::G4Cons - invalid pDPhi") ;
                                                   >> 121   
                                                   >> 122 // Ensure pSPhi in 0-2PI or -2PI-0 range if shape crosses 0
                                                   >> 123 
                                                   >> 124     if ( pSPhi < 0 ) fSPhi = 2.0*M_PI - fmod(fabs(pSPhi),2.0*M_PI) ;
                                                   >> 125     else             fSPhi = fmod(pSPhi,2.0*M_PI) ;
                                                   >> 126       
                                                   >> 127     if (fSPhi + fDPhi > 2.0*M_PI) fSPhi -= 2.0*M_PI ;
                                                   >> 128   }
175 }                                                 129 }
176                                                   130 
177 //////////////////////////////////////////////    131 /////////////////////////////////////////////////////////////////////
178 //                                                132 //
179 // Return whether point inside/outside/on surf    133 // Return whether point inside/outside/on surface
180                                                   134 
181 EInside G4Cons::Inside(const G4ThreeVector& p)    135 EInside G4Cons::Inside(const G4ThreeVector& p) const
182 {                                                 136 {
183   G4double r2, rl, rh, pPhi, tolRMin, tolRMax; << 137   G4double r2, rl, rh, pPhi, tolRMin, tolRMax ;
184   EInside in;                                  << 138   EInside in = kOutside ;
185                                                << 
186   if (std::fabs(p.z()) > fDz + halfCarToleranc << 
187   else if(std::fabs(p.z()) >= fDz - halfCarTol << 
188   else                                         << 
189                                                   139 
190   r2 = p.x()*p.x() + p.y()*p.y() ;             << 140   if (fabs(p.z()) <= fDz - kCarTolerance*0.5 )
191   rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz << 141   {
192   rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z << 142     r2 = p.x()*p.x() + p.y()*p.y() ;
                                                   >> 143     rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
193                                                   144 
194   // rh2 = rh*rh;                              << 145 // inner radius at z of point
195                                                   146 
196   tolRMin = rl - halfRadTolerance;             << 147     rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
197   if ( tolRMin < 0 )  { tolRMin = 0; }         << 148 // outer radius at z of point
198   tolRMax = rh + halfRadTolerance;             << 
199                                                   149 
200   if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tol << 
201                                                   150 
202   if (rl != 0.0) { tolRMin = rl + halfRadToler << 151     if (rl) tolRMin = rl + kRadTolerance*0.5 ;
203   else    { tolRMin = 0.0; }                   << 152     else    tolRMin = 0.0 ;
204   tolRMax = rh - halfRadTolerance;             << 
205                                                << 
206   if (in == kInside) // else it's kSurface alr << 
207   {                                            << 
208      if ( (r2 < tolRMin*tolRMin) || (r2 >= tol << 
209   }                                            << 
210   if ( !fPhiFullCone && ((p.x() != 0.0) || (p. << 
211   {                                            << 
212     pPhi = std::atan2(p.y(),p.x()) ;           << 
213                                                   153 
214     if ( pPhi < fSPhi - halfAngTolerance  )    << 154     tolRMax = rh - kRadTolerance*0.5 ;
215     else if ( pPhi > fSPhi + fDPhi + halfAngTo << 155       
216                                                << 156     if ( r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax )
217     if ( (pPhi < fSPhi - halfAngTolerance) ||  << 
218          (pPhi > fSPhi + fDPhi + halfAngTolera << 
219                                                << 
220     else if (in == kInside)  // else it's kSur << 
221     {                                             157     {
222        if ( (pPhi < fSPhi + halfAngTolerance)  << 158       if ( fDPhi == 2*M_PI || r2 == 0 ) in = kInside ;
223             (pPhi > fSPhi + fDPhi - halfAngTol << 159       else
                                                   >> 160       {
                                                   >> 161 // Try inner tolerant phi boundaries (=>inside)
                                                   >> 162 // if not inside, try outer tolerant phi boundaries
                                                   >> 163 
                                                   >> 164   pPhi = atan2(p.y(),p.x()) ;
                                                   >> 165 
                                                   >> 166   if (pPhi < 0 )   pPhi += 2*M_PI ; // 0<=pPhi<2pi
                                                   >> 167   if (fSPhi >= 0)
                                                   >> 168   {
                                                   >> 169     if ( pPhi >= fSPhi + kAngTolerance*0.5 &&
                                                   >> 170          pPhi <= fSPhi + fDPhi - kAngTolerance*0.5  )
                                                   >> 171     {
                                                   >> 172       in = kInside ;
                                                   >> 173     }
                                                   >> 174     else if ( pPhi >= fSPhi - kAngTolerance*0.5 &&
                                                   >> 175         pPhi <= fSPhi + fDPhi + kAngTolerance*0.5 )
                                                   >> 176     {
                                                   >> 177       in = kSurface ;
                                                   >> 178     }
                                                   >> 179   }
                                                   >> 180   else
                                                   >> 181   {
                                                   >> 182     if ( pPhi < fSPhi + 2*M_PI) pPhi += 2*M_PI ;
                                                   >> 183 
                                                   >> 184     if ( pPhi >= fSPhi + 2*M_PI + kAngTolerance*0.5 &&
                                                   >> 185          pPhi <= fSPhi + fDPhi + 2*M_PI - kAngTolerance*0.5   )
                                                   >> 186     {
                                                   >> 187       in = kInside ;
                                                   >> 188     }
                                                   >> 189     else if ( pPhi >= fSPhi + 2*M_PI - kAngTolerance*0.5 &&
                                                   >> 190         pPhi <= fSPhi + fDPhi + 2*M_PI + kAngTolerance*0.5  )
                                                   >> 191     {
                                                   >> 192       in = kSurface ;
                                                   >> 193     }
                                                   >> 194   }       
                                                   >> 195       }
                                                   >> 196     }
                                                   >> 197     else   // Try generous boundaries
                                                   >> 198     {
                                                   >> 199       tolRMin = rl - kRadTolerance*0.5 ;
                                                   >> 200       tolRMax = rh + kRadTolerance*0.5 ;
                                                   >> 201 
                                                   >> 202       if ( tolRMin < 0 ) tolRMin = 0 ;
                                                   >> 203 
                                                   >> 204       if (r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
                                                   >> 205       {
                                                   >> 206   if (fDPhi == 2*M_PI || r2 == 0) // Continuous in phi or on z-axis
                                                   >> 207         {
                                                   >> 208     in = kSurface ;
                                                   >> 209   }
                                                   >> 210   else  // Try outer tolerant phi boundaries only
                                                   >> 211   {
                                                   >> 212     pPhi = atan2(p.y(),p.x()) ;
                                                   >> 213 
                                                   >> 214     if (pPhi < 0) pPhi += 2*M_PI ;  // 0<=pPhi<2pi
                                                   >> 215 
                                                   >> 216     if ( fSPhi >= 0 )
                                                   >> 217     {
                                                   >> 218       if ( pPhi >= fSPhi - kAngTolerance*0.5 &&
                                                   >> 219      pPhi <= fSPhi + fDPhi + kAngTolerance*0.5  )
                                                   >> 220       {
                                                   >> 221         in = kSurface ;
                                                   >> 222       }
                                                   >> 223     }
                                                   >> 224     else
                                                   >> 225     {
                                                   >> 226       if ( pPhi < fSPhi + 2*M_PI) pPhi += 2*M_PI ;
                                                   >> 227 
                                                   >> 228       if ( pPhi >= fSPhi + 2*M_PI - kAngTolerance*0.5 &&
                                                   >> 229      pPhi <= fSPhi + fDPhi + 2*M_PI + kAngTolerance*0.5  )
                                                   >> 230       {
                                                   >> 231         in = kSurface ;
                                                   >> 232       }
                                                   >> 233     }
                                                   >> 234   }
                                                   >> 235       }
                                                   >> 236     }
                                                   >> 237   }
                                                   >> 238   else if (fabs(p.z()) <= fDz + kCarTolerance*0.5 )
                                                   >> 239   {
                                                   >> 240 // Check within tolerant r limits
                                                   >> 241 
                                                   >> 242     r2 = p.x()*p.x() + p.y()*p.y() ;
                                                   >> 243     rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
                                                   >> 244 
                                                   >> 245 // inner radius at z of point
                                                   >> 246 
                                                   >> 247     rh = 0.5*(fRmax2*(p.z() + fDz) + fRmax1*(fDz - p.z()))/fDz;
                                                   >> 248 
                                                   >> 249 // outer radius at z of point
                                                   >> 250 
                                                   >> 251     tolRMin = rl - kRadTolerance*0.5 ;
                                                   >> 252 
                                                   >> 253     if (tolRMin < 0 ) tolRMin = 0 ;
                                                   >> 254 
                                                   >> 255     tolRMax = rh + kRadTolerance*0.5 ;
                                                   >> 256 
                                                   >> 257     if (r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
                                                   >> 258     {
                                                   >> 259       if ( fDPhi == 2*M_PI || r2 == 0 ) // Continuous in phi or on z-axis
                                                   >> 260       {
                                                   >> 261   in = kSurface ;
                                                   >> 262       }
                                                   >> 263       else // Try outer tolerant phi boundaries
                                                   >> 264       {
                                                   >> 265   pPhi = atan2(p.y(),p.x()) ;
                                                   >> 266 
                                                   >> 267   if ( pPhi  <  0 )  pPhi += 2*M_PI ;   // 0<=pPhi<2pi
                                                   >> 268   if ( fSPhi >= 0 )
                                                   >> 269   {
                                                   >> 270     if ( pPhi >= fSPhi - kAngTolerance*0.5 &&
                                                   >> 271          pPhi <= fSPhi + fDPhi + kAngTolerance*0.5 )
                                                   >> 272     {
                                                   >> 273       in = kSurface ;
                                                   >> 274     }
                                                   >> 275   }
                                                   >> 276   else
                                                   >> 277   {
                                                   >> 278     if ( pPhi <  fSPhi + 2*M_PI) pPhi += 2*M_PI ;
                                                   >> 279     if ( pPhi >= fSPhi + 2*M_PI - kAngTolerance*0.5 &&
                                                   >> 280          pPhi <= fSPhi + fDPhi + 2*M_PI + kAngTolerance*0.5 )
                                                   >> 281     {
                                                   >> 282       in = kSurface ;
                                                   >> 283     }
                                                   >> 284   }    
                                                   >> 285       }
224     }                                             286     }
225   }                                               287   }
226   else if ( !fPhiFullCone )  { in = kSurface;  << 
227                                                << 
228   return in ;                                     288   return in ;
229 }                                                 289 }
230                                                   290 
231 //////////////////////////////////////////////    291 /////////////////////////////////////////////////////////////////////////
232 //                                                292 //
233 // Dispatch to parameterisation for replicatio    293 // Dispatch to parameterisation for replication mechanism dimension
234 // computation & modification.                    294 // computation & modification.
235                                                   295 
236 void G4Cons::ComputeDimensions(      G4VPVPara    296 void G4Cons::ComputeDimensions(      G4VPVParameterisation* p,
237                                const G4int        297                                const G4int                  n,
238                                const G4VPhysic    298                                const G4VPhysicalVolume*     pRep    )
239 {                                                 299 {
240   p->ComputeDimensions(*this,n,pRep) ;            300   p->ComputeDimensions(*this,n,pRep) ;
241 }                                                 301 }
242                                                   302 
243 ////////////////////////////////////////////// << 
244 //                                             << 
245 // Get bounding box                            << 
246                                                   303 
247 void G4Cons::BoundingLimits(G4ThreeVector& pMi << 304 ///////////////////////////////////////////////////////////////////////////
248 {                                              << 
249   G4double rmin = std::min(GetInnerRadiusMinus << 
250   G4double rmax = std::max(GetOuterRadiusMinus << 
251   G4double dz   = GetZHalfLength();            << 
252                                                << 
253   // Find bounding box                         << 
254   //                                           << 
255   if (GetDeltaPhiAngle() < twopi)              << 
256   {                                            << 
257     G4TwoVector vmin,vmax;                     << 
258     G4GeomTools::DiskExtent(rmin,rmax,         << 
259                             GetSinStartPhi(),G << 
260                             GetSinEndPhi(),Get << 
261                             vmin,vmax);        << 
262     pMin.set(vmin.x(),vmin.y(),-dz);           << 
263     pMax.set(vmax.x(),vmax.y(), dz);           << 
264   }                                            << 
265   else                                         << 
266   {                                            << 
267     pMin.set(-rmax,-rmax,-dz);                 << 
268     pMax.set( rmax, rmax, dz);                 << 
269   }                                            << 
270                                                << 
271   // Check correctness of the bounding box     << 
272   //                                           << 
273   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
274   {                                            << 
275     std::ostringstream message;                << 
276     message << "Bad bounding box (min >= max)  << 
277             << GetName() << " !"               << 
278             << "\npMin = " << pMin             << 
279             << "\npMax = " << pMax;            << 
280     G4Exception("G4Cons::BoundingLimits()", "G << 
281                 JustWarning, message);         << 
282     DumpInfo();                                << 
283   }                                            << 
284 }                                              << 
285                                                << 
286 ////////////////////////////////////////////// << 
287 //                                                305 //
288 // Calculate extent under transform and specif    306 // Calculate extent under transform and specified limit
289                                                   307 
290 G4bool G4Cons::CalculateExtent( const EAxis       308 G4bool G4Cons::CalculateExtent( const EAxis              pAxis,
291                                 const G4VoxelL << 309               const G4VoxelLimits&     pVoxelLimit,
292                                 const G4Affine << 310               const G4AffineTransform& pTransform,
293                                       G4double << 311                     G4double&          pMin,
294                                       G4double << 312                     G4double&          pMax  ) const
295 {                                                 313 {
296   G4ThreeVector bmin, bmax;                    << 314   if ( !pTransform.IsRotated() && 
297   G4bool exist;                                << 315         fDPhi == 2.0*M_PI && fRmin1 == 0 && fRmin2 == 0 )
298                                                << 
299   // Get bounding box                          << 
300   BoundingLimits(bmin,bmax);                   << 
301                                                << 
302   // Check bounding box                        << 
303   G4BoundingEnvelope bbox(bmin,bmax);          << 
304 #ifdef G4BBOX_EXTENT                           << 
305   if (true) return bbox.CalculateExtent(pAxis, << 
306 #endif                                         << 
307   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 
308   {                                               316   {
309     return exist = pMin < pMax;                << 317 // Special case handling for unrotated solid cones
310   }                                            << 318 // Compute z/x/y mins and maxs for bounding box respecting limits,
311                                                << 319 // with early returns if outside limits. Then switch() on pAxis,
312   // Get parameters of the solid               << 320 // and compute exact x and y limit for x/y case
313   G4double rmin1 = GetInnerRadiusMinusZ();     << 321       
314   G4double rmax1 = GetOuterRadiusMinusZ();     << 322     G4double xoffset, xMin, xMax ;
315   G4double rmin2 = GetInnerRadiusPlusZ();      << 323     G4double yoffset, yMin, yMax ;
316   G4double rmax2 = GetOuterRadiusPlusZ();      << 324     G4double zoffset, zMin, zMax ;
317   G4double dz    = GetZHalfLength();           << 325 
318   G4double dphi  = GetDeltaPhiAngle();         << 326     G4double diff1, diff2, maxDiff, newMin, newMax, RMax ;
319                                                << 327     G4double xoff1, xoff2, yoff1, yoff2 ;
320   // Find bounding envelope and calculate exte << 328       
321   //                                           << 329     zoffset = pTransform.NetTranslation().z();
322   const G4int NSTEPS = 24;            // numbe << 330     zMin    = zoffset - fDz ;
323   G4double astep  = twopi/NSTEPS;     // max a << 331     zMax    = zoffset + fDz ;
324   G4int    ksteps = (dphi <= astep) ? 1 : (G4i << 
325   G4double ang    = dphi/ksteps;               << 
326                                                << 
327   G4double sinHalf = std::sin(0.5*ang);        << 
328   G4double cosHalf = std::cos(0.5*ang);        << 
329   G4double sinStep = 2.*sinHalf*cosHalf;       << 
330   G4double cosStep = 1. - 2.*sinHalf*sinHalf;  << 
331   G4double rext1   = rmax1/cosHalf;            << 
332   G4double rext2   = rmax2/cosHalf;            << 
333                                                << 
334   // bounding envelope for full cone without h << 
335   // in other cases it is a sequence of quadri << 
336   if (rmin1 == 0 && rmin2 == 0 && dphi == twop << 
337   {                                            << 
338     G4double sinCur = sinHalf;                 << 
339     G4double cosCur = cosHalf;                 << 
340                                                << 
341     G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 
342     for (G4int k=0; k<NSTEPS; ++k)             << 
343     {                                          << 
344       baseA[k].set(rext1*cosCur,rext1*sinCur,- << 
345       baseB[k].set(rext2*cosCur,rext2*sinCur,  << 
346                                                << 
347       G4double sinTmp = sinCur;                << 
348       sinCur = sinCur*cosStep + cosCur*sinStep << 
349       cosCur = cosCur*cosStep - sinTmp*sinStep << 
350     }                                          << 
351     std::vector<const G4ThreeVectorList *> pol << 
352     polygons[0] = &baseA;                      << 
353     polygons[1] = &baseB;                      << 
354     G4BoundingEnvelope benv(bmin,bmax,polygons << 
355     exist = benv.CalculateExtent(pAxis,pVoxelL << 
356   }                                            << 
357   else                                         << 
358   {                                            << 
359     G4double sinStart = GetSinStartPhi();      << 
360     G4double cosStart = GetCosStartPhi();      << 
361     G4double sinEnd   = GetSinEndPhi();        << 
362     G4double cosEnd   = GetCosEndPhi();        << 
363     G4double sinCur   = sinStart*cosHalf + cos << 
364     G4double cosCur   = cosStart*cosHalf - sin << 
365                                                << 
366     // set quadrilaterals                      << 
367     G4ThreeVectorList pols[NSTEPS+2];          << 
368     for (G4int k=0; k<ksteps+2; ++k) pols[k].r << 
369     pols[0][0].set(rmin2*cosStart,rmin2*sinSta << 
370     pols[0][1].set(rmin1*cosStart,rmin1*sinSta << 
371     pols[0][2].set(rmax1*cosStart,rmax1*sinSta << 
372     pols[0][3].set(rmax2*cosStart,rmax2*sinSta << 
373     for (G4int k=1; k<ksteps+1; ++k)           << 
374     {                                          << 
375       pols[k][0].set(rmin2*cosCur,rmin2*sinCur << 
376       pols[k][1].set(rmin1*cosCur,rmin1*sinCur << 
377       pols[k][2].set(rext1*cosCur,rext1*sinCur << 
378       pols[k][3].set(rext2*cosCur,rext2*sinCur << 
379                                                << 
380       G4double sinTmp = sinCur;                << 
381       sinCur = sinCur*cosStep + cosCur*sinStep << 
382       cosCur = cosCur*cosStep - sinTmp*sinStep << 
383     }                                          << 
384     pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*s << 
385     pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*s << 
386     pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*s << 
387     pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*s << 
388                                                << 
389     // set envelope and calculate extent       << 
390     std::vector<const G4ThreeVectorList *> pol << 
391     polygons.resize(ksteps+2);                 << 
392     for (G4int k=0; k<ksteps+2; ++k) polygons[ << 
393     G4BoundingEnvelope benv(bmin,bmax,polygons << 
394     exist = benv.CalculateExtent(pAxis,pVoxelL << 
395   }                                            << 
396   return exist;                                << 
397 }                                              << 
398                                                << 
399 ////////////////////////////////////////////// << 
400 //                                             << 
401 // Return unit normal of surface closest to p  << 
402 // - note if point on z axis, ignore phi divid << 
403 // - unsafe if point close to z axis a rmin=0  << 
404                                                   332 
405 G4ThreeVector G4Cons::SurfaceNormal( const G4T << 333     if (pVoxelLimit.IsZLimited())
406 {                                              << 334     {
407   G4int noSurfaces = 0;                        << 335       if( zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance || 
408   G4double rho, pPhi;                          << 336           zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance   )
409   G4double distZ, distRMin, distRMax;          << 337       {
410   G4double distSPhi = kInfinity, distEPhi = kI << 338   return false ;
411   G4double tanRMin, secRMin, pRMin, widRMin;   << 339       }
412   G4double tanRMax, secRMax, pRMax, widRMax;   << 340       else
                                                   >> 341       {
                                                   >> 342   if ( zMin < pVoxelLimit.GetMinZExtent() )
                                                   >> 343   {
                                                   >> 344     zMin = pVoxelLimit.GetMinZExtent() ;
                                                   >> 345   }
                                                   >> 346   if ( zMax > pVoxelLimit.GetMaxZExtent() )
                                                   >> 347   {
                                                   >> 348     zMax = pVoxelLimit.GetMaxZExtent() ;
                                                   >> 349   }
                                                   >> 350       }
                                                   >> 351     }
                                                   >> 352     xoffset = pTransform.NetTranslation().x() ;
                                                   >> 353     RMax    = (fRmax2 >= fRmax1) ?  zMax : zMin  ;                  
                                                   >> 354     xMax    = xoffset + (fRmax1 + fRmax2)*0.5 + 
                                                   >> 355               (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
                                                   >> 356     xMin    = 2*xoffset-xMax ;
413                                                   357 
414   G4ThreeVector norm, sumnorm(0.,0.,0.), nZ =  << 358     if (pVoxelLimit.IsXLimited())
415   G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;    << 359     {
                                                   >> 360       if ( xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance || 
                                                   >> 361            xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance    )
                                                   >> 362       {
                                                   >> 363   return false ;
                                                   >> 364       }
                                                   >> 365       else
                                                   >> 366       {
                                                   >> 367   if ( xMin < pVoxelLimit.GetMinXExtent() )
                                                   >> 368   {
                                                   >> 369     xMin = pVoxelLimit.GetMinXExtent() ;
                                                   >> 370   }
                                                   >> 371   if ( xMax > pVoxelLimit.GetMaxXExtent() )
                                                   >> 372   {
                                                   >> 373             xMax=pVoxelLimit.GetMaxXExtent() ;
                                                   >> 374   }
                                                   >> 375       }
                                                   >> 376     }
                                                   >> 377     yoffset = pTransform.NetTranslation().y() ;
                                                   >> 378     yMax    = yoffset + (fRmax1 + fRmax2)*0.5 + 
                                                   >> 379               (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
                                                   >> 380     yMin    = 2*yoffset-yMax ;
                                                   >> 381     RMax    = yMax - yoffset ;  // = max radius due to Zmax/Zmin cuttings
416                                                   382 
417   distZ = std::fabs(std::fabs(p.z()) - fDz);   << 383     if (pVoxelLimit.IsYLimited())
418   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) << 384     {
                                                   >> 385       if ( yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance || 
                                                   >> 386            yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance     )
                                                   >> 387       {
                                                   >> 388   return false ;
                                                   >> 389       }
                                                   >> 390       else
                                                   >> 391       {
                                                   >> 392   if ( yMin < pVoxelLimit.GetMinYExtent() )
                                                   >> 393   {
                                                   >> 394     yMin = pVoxelLimit.GetMinYExtent() ;
                                                   >> 395   }
                                                   >> 396   if ( yMax > pVoxelLimit.GetMaxYExtent() )
                                                   >> 397   {
                                                   >> 398     yMax = pVoxelLimit.GetMaxYExtent() ;
                                                   >> 399   }
                                                   >> 400       }
                                                   >> 401     }    
                                                   >> 402     switch (pAxis) // Known to cut cones
                                                   >> 403     {
                                                   >> 404       case kXAxis:
                                                   >> 405   yoff1 = yoffset - yMin ;
                                                   >> 406   yoff2 = yMax - yoffset ;
                                                   >> 407 
                                                   >> 408   if (yoff1 >= 0 && yoff2 >= 0) // Y limits cross max/min x => no change
                                                   >> 409   {
                                                   >> 410     pMin = xMin ;
                                                   >> 411     pMax = xMax ;
                                                   >> 412   }
                                                   >> 413   else
                                                   >> 414   {
                                                   >> 415 // Y limits don't cross max/min x => compute max delta x, hence new mins/maxs
                                                   >> 416          
                                                   >> 417     diff1   = sqrt(RMax*RMax - yoff1*yoff1) ;
                                                   >> 418     diff2   = sqrt(RMax*RMax - yoff2*yoff2) ;
                                                   >> 419     maxDiff = (diff1>diff2) ? diff1:diff2 ;
                                                   >> 420     newMin  = xoffset - maxDiff ;
                                                   >> 421     newMax  = xoffset + maxDiff ;
                                                   >> 422     pMin    = ( newMin < xMin ) ? xMin : newMin  ;
                                                   >> 423     pMax    = ( newMax > xMax) ? xMax : newMax ;
                                                   >> 424   } 
                                                   >> 425       break ;
419                                                   426 
420   tanRMin  = (fRmin2 - fRmin1)*0.5/fDz;        << 427       case kYAxis:
421   secRMin  = std::sqrt(1 + tanRMin*tanRMin);   << 428   xoff1 = xoffset - xMin ;
422   pRMin    = rho - p.z()*tanRMin;              << 429   xoff2 = xMax - xoffset ;
423   widRMin  = fRmin2 - fDz*tanRMin;             << 430 
424   distRMin = std::fabs(pRMin - widRMin)/secRMi << 431   if (xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change
                                                   >> 432   {
                                                   >> 433     pMin = yMin ;
                                                   >> 434     pMax = yMax ;
                                                   >> 435   }
                                                   >> 436   else
                                                   >> 437   {
                                                   >> 438 // X limits don't cross max/min y => compute max delta y, hence new mins/maxs
                                                   >> 439 
                                                   >> 440     diff1   = sqrt(RMax*RMax - xoff1*xoff1) ;
                                                   >> 441     diff2   = sqrt(RMax*RMax-xoff2*xoff2) ;
                                                   >> 442     maxDiff = (diff1 > diff2) ? diff1:diff2 ;
                                                   >> 443     newMin  = yoffset - maxDiff ;
                                                   >> 444     newMax  = yoffset + maxDiff ;
                                                   >> 445     pMin    = (newMin < yMin) ? yMin : newMin ;
                                                   >> 446     pMax    = (newMax > yMax) ? yMax : newMax ;
                                                   >> 447   }
                                                   >> 448       break ;
425                                                   449 
426   tanRMax  = (fRmax2 - fRmax1)*0.5/fDz;        << 450       case kZAxis:
427   secRMax  = std::sqrt(1+tanRMax*tanRMax);     << 451   pMin = zMin ;
428   pRMax    = rho - p.z()*tanRMax;              << 452   pMax = zMax ;
429   widRMax  = fRmax2 - fDz*tanRMax;             << 453       break ;
430   distRMax = std::fabs(pRMax - widRMax)/secRMa << 454       
                                                   >> 455       default:
                                                   >> 456       break ;
                                                   >> 457     }
                                                   >> 458     pMin -= kCarTolerance ;
                                                   >> 459     pMax += kCarTolerance ;
431                                                   460 
432   if (!fPhiFullCone)   // Protected against (0 << 461     return true ;
                                                   >> 462   }
                                                   >> 463   else   // Calculate rotated vertex coordinates
433   {                                               464   {
434     if ( rho != 0.0 )                          << 465     G4int i, noEntries, noBetweenSections4 ;
435     {                                          << 466     G4bool existsAfterClip = false ;
436       pPhi = std::atan2(p.y(),p.x());          << 467     G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
437                                                   468 
438       if (pPhi  < fSPhi-halfCarTolerance)      << 469     pMin = +kInfinity ;
439       else if (pPhi > fSPhi+fDPhi+halfCarToler << 470     pMax = -kInfinity ;
440                                                   471 
441       distSPhi = std::fabs( pPhi - fSPhi );    << 472     noEntries          = vertices->size() ;
442       distEPhi = std::fabs( pPhi - fSPhi - fDP << 473     noBetweenSections4 = noEntries-4 ;
443     }                                          << 474       
444     else if( ((fRmin1) == 0.0) || ((fRmin2) == << 475     for ( i = 0 ; i < noEntries ; i += 4 )
445     {                                             476     {
446       distSPhi = 0.;                           << 477       ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
447       distEPhi = 0.;                           << 
448     }                                             478     }
449     nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0  << 479     for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
450     nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0  << 
451   }                                            << 
452   if ( rho > halfCarTolerance )                << 
453   {                                            << 
454     nR = G4ThreeVector(p.x()/rho/secRMax, p.y( << 
455     if ((fRmin1 != 0.0) || (fRmin2 != 0.0))    << 
456     {                                             480     {
457       nr = G4ThreeVector(-p.x()/rho/secRMin,-p << 481       ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
458     }                                          << 482     }    
459   }                                            << 483     if ( pMin != kInfinity || pMax != -kInfinity )
460                                                << 
461   if( distRMax <= halfCarTolerance )           << 
462   {                                            << 
463     ++noSurfaces;                              << 
464     sumnorm += nR;                             << 
465   }                                            << 
466   if( ((fRmin1 != 0.0) || (fRmin2 != 0.0)) &&  << 
467   {                                            << 
468     ++noSurfaces;                              << 
469     sumnorm += nr;                             << 
470   }                                            << 
471   if( !fPhiFullCone )                          << 
472   {                                            << 
473     if (distSPhi <= halfAngTolerance)          << 
474     {                                             484     {
475       ++noSurfaces;                            << 485       existsAfterClip = true ;
476       sumnorm += nPs;                          << 486         
                                                   >> 487 // Add 2*tolerance to avoid precision troubles
                                                   >> 488 
                                                   >> 489       pMin -= kCarTolerance ;
                                                   >> 490       pMax += kCarTolerance ;
477     }                                             491     }
478     if (distEPhi <= halfAngTolerance)          << 492     else
479     {                                             493     {
480       ++noSurfaces;                            << 494 // Check for case where completely enveloping clipping volume
481       sumnorm += nPe;                          << 495 // If point inside then we are confident that the solid completely
                                                   >> 496 // envelopes the clipping volume. Hence set min/max extents according
                                                   >> 497 // to clipping volume extents along the specified axis.
                                                   >> 498        
                                                   >> 499       G4ThreeVector clipCentre(
                                                   >> 500       (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
                                                   >> 501       (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
                                                   >> 502       (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5  ) ;
                                                   >> 503         
                                                   >> 504       if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
                                                   >> 505       {
                                                   >> 506   existsAfterClip = true ;
                                                   >> 507   pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
                                                   >> 508         pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
                                                   >> 509       }
482     }                                             510     }
                                                   >> 511     delete vertices ;
                                                   >> 512     return existsAfterClip ;
483   }                                               513   }
484   if (distZ <= halfCarTolerance)               << 
485   {                                            << 
486     ++noSurfaces;                              << 
487     if ( p.z() >= 0.)  { sumnorm += nZ; }      << 
488     else               { sumnorm -= nZ; }      << 
489   }                                            << 
490   if ( noSurfaces == 0 )                       << 
491   {                                            << 
492 #ifdef G4CSGDEBUG                              << 
493     G4Exception("G4Cons::SurfaceNormal(p)", "G << 
494                 JustWarning, "Point p is not o << 
495 #endif                                         << 
496      norm = ApproxSurfaceNormal(p);            << 
497   }                                            << 
498   else if ( noSurfaces == 1 )  { norm = sumnor << 
499   else                         { norm = sumnor << 
500                                                << 
501   return norm ;                                << 
502 }                                                 514 }
503                                                   515 
504 ////////////////////////////////////////////// << 516 ////////////////////////////////////////////////////////////////////////
505 //                                                517 //
506 // Algorithm for SurfaceNormal() following the << 518 // Return unit normal of surface closest to p
507 // for points not on the surface               << 519 // - note if point on z axis, ignore phi divided sides
                                                   >> 520 // - unsafe if point close to z axis a rmin=0 - no explicit checks
508                                                   521 
509 G4ThreeVector G4Cons::ApproxSurfaceNormal( con << 522 G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const
510 {                                                 523 {
511   ENorm side ;                                    524   ENorm side ;
512   G4ThreeVector norm ;                            525   G4ThreeVector norm ;
513   G4double rho, phi ;                             526   G4double rho, phi ;
514   G4double distZ, distRMin, distRMax, distSPhi    527   G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
515   G4double tanRMin, secRMin, pRMin, widRMin ;     528   G4double tanRMin, secRMin, pRMin, widRMin ;
516   G4double tanRMax, secRMax, pRMax, widRMax ;     529   G4double tanRMax, secRMax, pRMax, widRMax ;
517                                                   530 
518   distZ = std::fabs(std::fabs(p.z()) - fDz) ;  << 531   distZ = fabs(fabs(p.z()) - fDz) ;
519   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) << 532   rho   = sqrt(p.x()*p.x() + p.y()*p.y()) ;
520                                                   533 
521   tanRMin  = (fRmin2 - fRmin1)*0.5/fDz ;          534   tanRMin  = (fRmin2 - fRmin1)*0.5/fDz ;
522   secRMin  = std::sqrt(1 + tanRMin*tanRMin) ;  << 535   secRMin  = sqrt(1 + tanRMin*tanRMin) ;
523   pRMin    = rho - p.z()*tanRMin ;                536   pRMin    = rho - p.z()*tanRMin ;
524   widRMin  = fRmin2 - fDz*tanRMin ;               537   widRMin  = fRmin2 - fDz*tanRMin ;
525   distRMin = std::fabs(pRMin - widRMin)/secRMi << 538   distRMin = fabs(pRMin - widRMin)/secRMin ;
526                                                   539 
527   tanRMax  = (fRmax2 - fRmax1)*0.5/fDz ;          540   tanRMax  = (fRmax2 - fRmax1)*0.5/fDz ;
528   secRMax  = std::sqrt(1+tanRMax*tanRMax) ;    << 541   secRMax  = sqrt(1+tanRMax*tanRMax) ;
529   pRMax    = rho - p.z()*tanRMax ;                542   pRMax    = rho - p.z()*tanRMax ;
530   widRMax  = fRmax2 - fDz*tanRMax ;               543   widRMax  = fRmax2 - fDz*tanRMax ;
531   distRMax = std::fabs(pRMax - widRMax)/secRMa << 544   distRMax = fabs(pRMax - widRMax)/secRMax ;
532                                                   545   
533   if (distRMin < distRMax)  // First minimum      546   if (distRMin < distRMax)  // First minimum
534   {                                               547   {
535     if (distZ < distRMin)                         548     if (distZ < distRMin)
536     {                                             549     {
537       distMin = distZ ;                           550       distMin = distZ ;
538       side    = kNZ ;                             551       side    = kNZ ;
539     }                                             552     }
540     else                                          553     else
541     {                                             554     {
542       distMin = distRMin ;                        555       distMin = distRMin ;
543       side    = kNRMin ;                          556       side    = kNRMin ;
544     }                                             557     }
545   }                                               558   }
546   else                                            559   else
547   {                                               560   {
548     if (distZ < distRMax)                         561     if (distZ < distRMax)
549     {                                             562     {
550       distMin = distZ ;                           563       distMin = distZ ;
551       side    = kNZ ;                             564       side    = kNZ ;
552     }                                             565     }
553     else                                          566     else
554     {                                             567     {
555       distMin = distRMax ;                        568       distMin = distRMax ;
556       side    = kNRMax ;                          569       side    = kNRMax ;
557     }                                             570     }
558   }                                               571   }
559   if ( !fPhiFullCone && (rho != 0.0) )  // Pro << 572   if ( fDPhi < 2.0*M_PI && rho )  // Protected against (0,0,z) 
560   {                                               573   {
561     phi = std::atan2(p.y(),p.x()) ;            << 574     phi = atan2(p.y(),p.x()) ;
562                                                   575 
563     if (phi < 0)  { phi += twopi; }            << 576     if (phi < 0) phi += 2*M_PI ;
564                                                   577 
565     if (fSPhi < 0)  { distSPhi = std::fabs(phi << 578     if (fSPhi < 0) distSPhi = fabs(phi - (fSPhi + 2.0*M_PI))*rho ;
566     else            { distSPhi = std::fabs(phi << 579     else           distSPhi = fabs(phi - fSPhi)*rho ;
567                                                   580 
568     distEPhi = std::fabs(phi - fSPhi - fDPhi)* << 581     distEPhi = fabs(phi - fSPhi - fDPhi)*rho ;
569                                                   582 
570     // Find new minimum                        << 583 // Find new minimum
571                                                   584 
572     if (distSPhi < distEPhi)                      585     if (distSPhi < distEPhi)
573     {                                             586     {
574       if (distSPhi < distMin)  { side = kNSPhi << 587       if (distSPhi < distMin) side = kNSPhi ;
575     }                                             588     }
576     else                                          589     else 
577     {                                             590     {
578       if (distEPhi < distMin)  { side = kNEPhi << 591       if (distEPhi < distMin) side = kNEPhi ;
579     }                                             592     }
580   }                                            << 593   }   
581   switch (side)                                   594   switch (side)
582   {                                               595   {
583     case kNRMin:      // Inner radius          << 596     case kNRMin:      // Inner radius
584     {                                          << 
585       rho *= secRMin ;                            597       rho *= secRMin ;
586       norm = G4ThreeVector(-p.x()/rho, -p.y()/ << 598       norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,tanRMin/secRMin) ;
587       break ;                                     599       break ;
588     }                                          << 600     case kNRMax:      // Outer radius
589     case kNRMax:      // Outer radius          << 
590     {                                          << 
591       rho *= secRMax ;                            601       rho *= secRMax ;
592       norm = G4ThreeVector(p.x()/rho, p.y()/rh << 602       norm = G4ThreeVector(p.x()/rho,p.y()/rho,-tanRMax/secRMax) ;
593       break ;                                     603       break ;
594     }                                          << 604     case kNZ:     // +/- dz
595     case kNZ:         // +/- dz                << 605       if (p.z() > 0) norm = G4ThreeVector(0,0,1) ; 
596     {                                          << 606       else           norm = G4ThreeVector(0,0,-1) ; 
597       if (p.z() > 0)  { norm = G4ThreeVector(0 << 
598       else            { norm = G4ThreeVector(0 << 
599       break ;                                     607       break ;
600     }                                          << 
601     case kNSPhi:                                  608     case kNSPhi:
602     {                                          << 609       norm = G4ThreeVector(sin(fSPhi),-cos(fSPhi),0) ;
603       norm = G4ThreeVector(sinSPhi, -cosSPhi,  << 
604       break ;                                     610       break ;
605     }                                          << 
606     case kNEPhi:                                  611     case kNEPhi:
607     {                                          << 612       norm=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0) ;
608       norm = G4ThreeVector(-sinEPhi, cosEPhi,  << 
609       break ;                                     613       break ;
610     }                                          << 614     default:
611     default:          // Should never reach th << 615       G4Exception("Logic error in G4Cons::SurfaceNormal") ;
612     {                                          << 
613       DumpInfo();                              << 
614       G4Exception("G4Cons::ApproxSurfaceNormal << 
615                   "GeomSolids1002", JustWarnin << 
616                   "Undefined side for valid su << 
617       break ;                                     616       break ;    
618     }                                          << 
619   }                                               617   }
620   return norm ;                                   618   return norm ;
621 }                                                 619 }
622                                                   620 
623 //////////////////////////////////////////////    621 ////////////////////////////////////////////////////////////////////////
624 //                                                622 //
625 // Calculate distance to shape from outside, a    623 // Calculate distance to shape from outside, along normalised vector
626 // - return kInfinity if no intersection, or i    624 // - return kInfinity if no intersection, or intersection distance <= tolerance
627 //                                                625 //
628 // - Compute the intersection with the z plane    626 // - Compute the intersection with the z planes 
629 //        - if at valid r, phi, return            627 //        - if at valid r, phi, return
630 //                                                628 //
631 // -> If point is outside cone, compute inters    629 // -> If point is outside cone, compute intersection with rmax1*0.5
632 //        - if at valid phi,z return              630 //        - if at valid phi,z return
633 //        - if inside outer cone, handle case     631 //        - if inside outer cone, handle case when on tolerant outer cone
634 //          boundary and heading inwards(->0 t    632 //          boundary and heading inwards(->0 to in)
635 //                                                633 //
636 // -> Compute intersection with inner cone, ta    634 // -> Compute intersection with inner cone, taking largest +ve root
637 //        - if valid (in z,phi), save intersct    635 //        - if valid (in z,phi), save intersction
638 //                                                636 //
639 //    -> If phi segmented, compute intersectio    637 //    -> If phi segmented, compute intersections with phi half planes
640 //        - return smallest of valid phi inter    638 //        - return smallest of valid phi intersections and
641 //          inner radius intersection             639 //          inner radius intersection
642 //                                                640 //
643 // NOTE:                                          641 // NOTE:
                                                   >> 642 // - Precalculations for phi trigonometry are Done `just in time'
644 // - `if valid' implies tolerant checking of i    643 // - `if valid' implies tolerant checking of intersection points
645 // - z, phi intersection from Tubs                644 // - z, phi intersection from Tubs
646                                                   645 
647 G4double G4Cons::DistanceToIn( const G4ThreeVe    646 G4double G4Cons::DistanceToIn( const G4ThreeVector& p,
648                                const G4ThreeVe << 647              const G4ThreeVector& v   ) const
649 {                                                 648 {
650   G4double snxt = kInfinity ;      // snxt = d << 649   G4double snxt = kInfinity ;     // snxt = default return value
651   const G4double dRmax = 50*(fRmax1+fRmax2);// << 650 
                                                   >> 651   G4bool seg ;        // true if segmented in phi
                                                   >> 652   G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0. ;
                                                   >> 653           // half dphi + outer tolerance
                                                   >> 654   G4double cPhi,sinCPhi=0.,cosCPhi=0. ; // central phi
652                                                   655 
653   G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;  / << 656   G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones
654   G4double tanRMin,secRMin,rMinAv,rMinOAv ;    << 657   G4double tanRMin,secRMin,rMinAv,rMinIAv,rMinOAv ;
655   G4double rout,rin ;                             658   G4double rout,rin ;
656                                                   659 
657   G4double tolORMin,tolORMin2,tolIRMin,tolIRMi << 660   G4double tolORMin,tolORMin2 ; // `generous' radii squared
658   G4double tolORMax2,tolIRMax,tolIRMax2 ;      << 661   G4double tolORMax2 ;
659   G4double tolODz,tolIDz ;                        662   G4double tolODz,tolIDz ;
660                                                   663 
661   G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2, << 664   G4double Dist,s,xi,yi,zi,ri=0.,rho2,cosPsi ; // Intersection point variables
662                                                   665 
663   G4double t1,t2,t3,b,c,d ;    // Quadratic so << 666   G4double t1,t2,t3,b,c,d ;   // Quadratic solver variables 
664   G4double nt1,nt2,nt3 ;                          667   G4double nt1,nt2,nt3 ;
665   G4double Comp ;                                 668   G4double Comp ;
                                                   >> 669   G4double cosSPhi,sinSPhi ;    // Trig for phi start intersect
                                                   >> 670   G4double ePhi,cosEPhi,sinEPhi ; // for phi end intersect
666                                                   671 
667   G4ThreeVector Normal;                        << 672 //
                                                   >> 673 // Set phi divided flag and precalcs
                                                   >> 674 //
                                                   >> 675   if (fDPhi < 2.0*M_PI)
                                                   >> 676   {
                                                   >> 677     seg        = true ;
                                                   >> 678     hDPhi      = 0.5*fDPhi ;    // half delta phi
                                                   >> 679     cPhi       = fSPhi + hDPhi ; ;
                                                   >> 680     hDPhiOT    = hDPhi + 0.5*kAngTolerance ;  // outers tol' half delta phi 
                                                   >> 681     hDPhiIT    = hDPhi - 0.5*kAngTolerance ;
                                                   >> 682     sinCPhi    = sin(cPhi) ;
                                                   >> 683     cosCPhi    = cos(cPhi) ;
                                                   >> 684     cosHDPhiOT = cos(hDPhiOT) ;
                                                   >> 685     cosHDPhiIT = cos(hDPhiIT) ;
                                                   >> 686   }
                                                   >> 687   else     seg = false ;
668                                                   688 
669   // Cone Precalcs                             << 689 // Cone Precalcs
670                                                   690 
671   tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;           691   tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
672   secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; << 692   secRMin = sqrt(1.0 + tanRMin*tanRMin) ;
673   rMinAv  = (fRmin1 + fRmin2)*0.5 ;               693   rMinAv  = (fRmin1 + fRmin2)*0.5 ;
674                                                   694 
675   if (rMinAv > halfRadTolerance)               << 695   if (rMinAv > kRadTolerance*0.5)
676   {                                               696   {
677     rMinOAv = rMinAv - halfRadTolerance ;      << 697     rMinOAv = rMinAv - kRadTolerance*0.5 ;
                                                   >> 698     rMinIAv = rMinAv + kRadTolerance*0.5 ;
678   }                                               699   }
679   else                                            700   else
680   {                                               701   {
681     rMinOAv = 0.0 ;                               702     rMinOAv = 0.0 ;
                                                   >> 703     rMinIAv = 0.0 ;
682   }                                               704   }  
683   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;           705   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
684   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; << 706   secRMax = sqrt(1.0 + tanRMax*tanRMax) ;
685   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;               707   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
686   rMaxOAv = rMaxAv + halfRadTolerance ;        << 708   rMaxOAv = rMaxAv + kRadTolerance*0.5 ;
687                                                   709    
688   // Intersection with z-surfaces              << 710 // Intersection with Z surfaces
689                                                   711 
690   tolIDz = fDz - halfCarTolerance ;            << 712   tolIDz = fDz - kCarTolerance*0.5 ;
691   tolODz = fDz + halfCarTolerance ;            << 713   tolODz = fDz + kCarTolerance*0.5 ;
692                                                   714 
693   if (std::fabs(p.z()) >= tolIDz)              << 715   if (fabs(p.z()) >= tolIDz)
694   {                                               716   {
695     if ( p.z()*v.z() < 0 )    // at +Z going i << 717     if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
696     {                                             718     {
697       sd = (std::fabs(p.z()) - fDz)/std::fabs( << 719       s = (fabs(p.z()) - fDz)/fabs(v.z()) ; // Z intersect distance
698                                                   720 
699       if( sd < 0.0 )  { sd = 0.0; }            << 721       if( s < 0.0 ) s = 0.0 ;                  // negative dist -> zero
700                                                   722 
701       xi   = p.x() + sd*v.x() ;  // Intersecti << 723       xi   = p.x() + s*v.x() ;  // Intersection coords
702       yi   = p.y() + sd*v.y() ;                << 724       yi   = p.y() + s*v.y() ;
703       rhoi2 = xi*xi + yi*yi  ;                 << 725       rho2 = xi*xi + yi*yi   ;
704                                                   726 
705       // Check validity of intersection        << 727 // Check validity of intersection
706       // Calculate (outer) tolerant radi^2 at  << 728 // Calculate (outer) tolerant radi^2 at intersecion
707                                                   729 
708       if (v.z() > 0)                              730       if (v.z() > 0)
709       {                                           731       {
710         tolORMin  = fRmin1 - halfRadTolerance* << 732   tolORMin  = fRmin1 - 0.5*kRadTolerance*secRMin ;
711         tolIRMin  = fRmin1 + halfRadTolerance* << 733   tolORMax2 = (fRmax1 + 0.5*kRadTolerance*secRMax)*(fRmax1 + 0.5*kRadTolerance*secRMax) ;
712         tolIRMax  = fRmax1 - halfRadTolerance* << 
713         // tolORMax2 = (fRmax1 + halfRadTolera << 
714         //             (fRmax1 + halfRadTolera << 
715       }                                           734       }
716       else                                        735       else
717       {                                           736       {
718         tolORMin  = fRmin2 - halfRadTolerance* << 737   tolORMin  = fRmin2 - 0.5*kRadTolerance*secRMin ;
719         tolIRMin  = fRmin2 + halfRadTolerance* << 738   tolORMax2 = (fRmax2 + 0.5*kRadTolerance*secRMax)*(fRmax2 + 0.5*kRadTolerance*secRMax) ;
720         tolIRMax  = fRmax2 - halfRadTolerance* << 
721         // tolORMax2 = (fRmax2 + halfRadTolera << 
722         //             (fRmax2 + halfRadTolera << 
723       }                                        << 
724       if ( tolORMin > 0 )                      << 
725       {                                        << 
726         // tolORMin2 = tolORMin*tolORMin ;     << 
727         tolIRMin2 = tolIRMin*tolIRMin ;        << 
728       }                                           739       }
729       else                                     << 740       if ( tolORMin > 0 ) tolORMin2 = tolORMin*tolORMin ;
730       {                                        << 741       else                tolORMin2 = 0.0 ;
731         // tolORMin2 = 0.0 ;                   << 742 
732         tolIRMin2 = 0.0 ;                      << 743       if (tolORMin2 <= rho2 && rho2 <= tolORMax2)
733       }                                        << 
734       if ( tolIRMax > 0 )  { tolIRMax2 = tolIR << 
735       else                 { tolIRMax2 = 0.0;  << 
736                                                << 
737       if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= t << 
738       {                                           744       {
739         if ( !fPhiFullCone && (rhoi2 != 0.0) ) << 745   if ( seg && rho2 )
740         {                                      << 746   {
741           // Psi = angle made with central (av << 747 // Psi = angle made with central (average) phi of shape
742                                                   748 
743           cosPsi = (xi*cosCPhi + yi*sinCPhi)/s << 749     cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rho2) ;
744                                                   750 
745           if (cosPsi >= cosHDPhiIT)  { return  << 751     if (cosPsi >= cosHDPhiOT) return s ;
746         }                                      << 752   }
747         else                                   << 753   else return s ;
748         {                                      << 
749           return sd;                           << 
750         }                                      << 
751       }                                           754       }
752     }                                             755     }
753     else  // On/outside extent, and heading aw    756     else  // On/outside extent, and heading away  -> cannot intersect
754     {                                             757     {
755       return snxt ;                            << 758       return snxt ; 
756     }                                             759     }
757   }                                               760   }
758                                                   761     
759 // ----> Can not intersect z surfaces             762 // ----> Can not intersect z surfaces
760                                                   763 
761                                                   764 
762 // Intersection with outer cone (possible retu    765 // Intersection with outer cone (possible return) and
763 //                   inner cone (must also che    766 //                   inner cone (must also check phi)
764 //                                                767 //
765 // Intersection point (xi,yi,zi) on line x=p.x    768 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
766 //                                                769 //
767 // Intersects with x^2+y^2=(a*z+b)^2              770 // Intersects with x^2+y^2=(a*z+b)^2
768 //                                                771 //
769 // where a=tanRMax or tanRMin                     772 // where a=tanRMax or tanRMin
770 //       b=rMaxAv  or rMinAv                      773 //       b=rMaxAv  or rMinAv
771 //                                                774 //
772 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a    775 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
773 //     t1                        t2               776 //     t1                        t2                      t3  
774 //                                                777 //
775 //  \--------u-------/       \-----------v----    778 //  \--------u-------/       \-----------v----------/ \---------w--------/
776 //                                                779 //
777                                                   780 
778   t1   = 1.0 - v.z()*v.z() ;                      781   t1   = 1.0 - v.z()*v.z() ;
779   t2   = p.x()*v.x() + p.y()*v.y() ;              782   t2   = p.x()*v.x() + p.y()*v.y() ;
780   t3   = p.x()*p.x() + p.y()*p.y() ;              783   t3   = p.x()*p.x() + p.y()*p.y() ;
781   rin  = tanRMin*p.z() + rMinAv ;                 784   rin  = tanRMin*p.z() + rMinAv ;
782   rout = tanRMax*p.z() + rMaxAv ;                 785   rout = tanRMax*p.z() + rMaxAv ;
783                                                   786 
784   // Outer Cone Intersection                   << 787 // Outer Cone Intersection
785   // Must be outside/on outer cone for valid i << 788 // Must be outside/on outer cone for valid intersection
786                                                   789 
787   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;    790   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
788   nt2 = t2 - tanRMax*v.z()*rout ;                 791   nt2 = t2 - tanRMax*v.z()*rout ;
789   nt3 = t3 - rout*rout ;                          792   nt3 = t3 - rout*rout ;
790                                                   793 
791   if (std::fabs(nt1) > kRadTolerance)  // Equa << 794   if (fabs(nt1) > kRadTolerance)  // Equation quadratic => 2 roots
792   {                                               795   {
793     b = nt2/nt1;                               << 796     if ( nt3 > rout*kRadTolerance*secRMax || rout < 0 )
794     c = nt3/nt1;                               << 
795     d = b*b-c  ;                               << 
796     if ( (nt3 > rout*rout*kRadTolerance*kRadTo << 
797       || (rout < 0) )                          << 
798     {                                             797     {
799       // If outside real cone (should be rho-r << 798 // If outside real cone (should be rho-rout>kRadTolerance*0.5 NOT rho^2 etc)
800       // NOT rho^2 etc) saves a std::sqrt() at << 799 // saves a sqrt() at expense of accuracy
                                                   >> 800 
                                                   >> 801       b = nt2/nt1 ;
                                                   >> 802       c = nt3/nt1 ;
                                                   >> 803       d = b*b-c   ;
801                                                   804 
802       if (d >= 0)                                 805       if (d >= 0)
803       {                                           806       {
804                                                << 807           
805         if ((rout < 0) && (nt3 <= 0))          << 808         if (rout < 0 && nt3 <= 0 )
806         {                                         809         {
807           // Inside `shadow cone' with -ve rad << 810 // Inside `shadow cone' with -ve radius -> 2nd root could be on real cone
808           // -> 2nd root could be on real cone << 
809                                                   811 
810           if (b>0) { sd = c/(-b-std::sqrt(d)); << 812           s = -b + sqrt(d) ;
811           else     { sd = -b + std::sqrt(d);   << 813   }
812         }                                      << 814   else
813         else                                   << 815         {
814         {                                      << 816     if ( b <= 0 && c >= 0 ) // both >=0, try smaller root
815           if ((b <= 0) && (c >= 0)) // both >= << 817     {
816           {                                    << 818       s = -b - sqrt(d) ;
817             sd=c/(-b+std::sqrt(d));            << 819     }
818           }                                    << 
819           else                                    820           else
820           {                                    << 821     {
821             if ( c <= 0 ) // second >=0        << 822       if ( c <= 0 ) // second >=0
822             {                                  << 823       {
823               sd = -b + std::sqrt(d) ;         << 824                       s = -b + sqrt(d) ;
824               if((sd<0.0) && (sd>-halfRadToler << 825       }
825             }                                  << 
826             else  // both negative, travel awa    826             else  // both negative, travel away
827             {                                  << 827       {
828               return kInfinity ;               << 828                       return kInfinity ;
829             }                                  << 829       }
830           }                                    << 830     }
831         }                                      << 831   }
832         if ( sd >= 0 )  // If 'forwards'. Chec << 832   if ( s > 0 )  // If 'forwards'. Check z intersection
833         {                                      << 833   {
834           if ( sd>dRmax ) // Avoid rounding er << 834     zi = p.z() + s*v.z() ;
835           {               // 64 bits systems.  << 835 
836             G4double fTerm = sd-std::fmod(sd,d << 836     if (fabs(zi) <= tolODz)
837             sd = fTerm + DistanceToIn(p+fTerm* << 837     {
838           }                                    << 838 // Z ok. Check phi intersection if reqd
839           zi = p.z() + sd*v.z() ;              << 839 
840                                                << 840       if ( ! seg )  return s ;
841           if (std::fabs(zi) <= tolODz)         << 841       else
842           {                                    << 842       {
843             // Z ok. Check phi intersection if << 843         xi     = p.x() + s*v.x() ;
844                                                << 844         yi     = p.y() + s*v.y() ;
845             if ( fPhiFullCone )  { return sd;  << 845         ri     = rMaxAv + zi*tanRMax ;
846             else                               << 846         cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
847             {                                  << 847 
848               xi     = p.x() + sd*v.x() ;      << 848         if ( cosPsi >= cosHDPhiOT ) return s ;
849               yi     = p.y() + sd*v.y() ;      << 849       }
850               ri     = rMaxAv + zi*tanRMax ;   << 850     }             
851               cosPsi = (xi*cosCPhi + yi*sinCPh << 851   }                // end if (s>0)
852                                                << 
853               if ( cosPsi >= cosHDPhiIT )  { r << 
854             }                                  << 
855           }                                    << 
856         }                // end if (sd>0)      << 
857       }                                           852       }
858     }                                             853     }
859     else                                          854     else
860     {                                             855     {
861       // Inside outer cone                     << 856 // Inside outer cone
862       // check not inside, and heading through << 857 // check not inside, and heading through G4Cons (-> 0 to in)
863                                                   858 
864       if ( ( t3  > (rin + halfRadTolerance*sec << 859       if ( t3  > (rin + kRadTolerance*0.5*secRMin)*
865                    (rin + halfRadTolerance*sec << 860                  (rin + kRadTolerance*0.5*secRMin) && 
866         && (nt2 < 0) && (d >= 0) && (std::fabs << 861            nt2 < 0                                 && 
867       {                                        << 862            fabs(p.z()) <= tolIDz )
868         // Inside cones, delta r -ve, inside z << 863       {
869         // Point is on the Surface => check Di << 864 // Inside cones, delta r -ve, inside z extent
870                                                << 865 
871         xi     = p.x() ;                       << 866         if (seg)
872         yi     = p.y()  ;                      << 
873         risec  = std::sqrt(xi*xi + yi*yi)*secR << 
874         Normal = G4ThreeVector(xi/risec,yi/ris << 
875         if ( !fPhiFullCone )                   << 
876         {                                         867         {
877           cosPsi = (p.x()*cosCPhi + p.y()*sinC << 868           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/sqrt(t3) ;
878           if ( cosPsi >= cosHDPhiIT )          << 869 
879           {                                    << 870     if (cosPsi >= cosHDPhiIT) return 0.0 ;
880             if ( Normal.dot(v) <= 0 )  { retur << 871   }
881           }                                    << 872   else  return 0.0 ;
882         }                                      << 
883         else                                   << 
884         {                                      << 
885           if ( Normal.dot(v) <= 0 )  { return  << 
886         }                                      << 
887       }                                           873       }
888     }                                             874     }
889   }                                               875   }
890   else  //  Single root case                      876   else  //  Single root case 
891   {                                               877   {
892     if ( std::fabs(nt2) > kRadTolerance )      << 878     if ( fabs(nt2) > kRadTolerance )
893     {                                             879     {
894       sd = -0.5*nt3/nt2 ;                      << 880       s = -0.5*nt3/nt2 ;
895                                                   881 
896       if ( sd < 0 )  { return kInfinity; }   / << 882       if ( s < 0 )   return kInfinity ;   // travel away
897       else  // sd >= 0,  If 'forwards'. Check  << 883       else  // s >= 0,  If 'forwards'. Check z intersection
898       {                                           884       {
899         zi = p.z() + sd*v.z() ;                << 885   zi = p.z() + s*v.z() ;
900                                                << 
901         if ((std::fabs(zi) <= tolODz) && (nt2  << 
902         {                                      << 
903           // Z ok. Check phi intersection if r << 
904                                                   886 
905           if ( fPhiFullCone )  { return sd; }  << 887   if (fabs(zi) <= tolODz && nt2 < 0)
906           else                                 << 888   {
907           {                                    << 889 // Z ok. Check phi intersection if reqd
908             xi     = p.x() + sd*v.x() ;        << 890 
909             yi     = p.y() + sd*v.y() ;        << 891     if ( ! seg ) return s ;
910             ri     = rMaxAv + zi*tanRMax ;     << 892     else
911             cosPsi = (xi*cosCPhi + yi*sinCPhi) << 893     {
912                                                << 894       xi     = p.x() + s*v.x() ;
913             if (cosPsi >= cosHDPhiIT)  { retur << 895       yi     = p.y() + s*v.y() ;
914           }                                    << 896       ri     = rMaxAv + zi*tanRMax ;
915         }                                      << 897       cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
916       }                                        << 898 
                                                   >> 899       if (cosPsi >= cosHDPhiOT) return s ;
                                                   >> 900     }
                                                   >> 901   }             
                                                   >> 902       }                
917     }                                             903     }
918     else  //    travel || cone surface from it    904     else  //    travel || cone surface from its origin
919     {                                             905     {
920       sd = kInfinity ;                         << 906       s = kInfinity ;
921     }                                             907     }
922   }                                               908   }
923                                                   909 
924   // Inner Cone Intersection                   << 910 // Inner Cone Intersection
925   // o Space is divided into 3 areas:          << 911 // o Space is divided into 3 areas:
926   //   1) Radius greater than real inner cone  << 912 //   1) Radius greater than real inner cone & imaginary cone & outside
927   //      tolerance                            << 913 //      tolerance
928   //   2) Radius less than inner or imaginary  << 914 //   2) Radius less than inner or imaginary cone & outside tolarance
929   //   3) Within tolerance of real or imaginar << 915 //   3) Within tolerance of real or imaginary cones
930   //      - Extra checks needed for 3's inters << 916 //      - Extra checks needed for 3's intersections => lots of duplicated code
931   //        => lots of duplicated code         << 
932                                                   917 
933   if (rMinAv != 0.0)                           << 918   if (rMinAv)
934   {                                            << 919   {
935     nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z())    920     nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
936     nt2 = t2 - tanRMin*v.z()*rin ;                921     nt2 = t2 - tanRMin*v.z()*rin ;
937     nt3 = t3 - rin*rin ;                          922     nt3 = t3 - rin*rin ;
938                                                   923  
939     if ( nt1 != 0.0 )                          << 924     if ( nt1 )
940     {                                             925     {
941       if ( nt3 > rin*kRadTolerance*secRMin )      926       if ( nt3 > rin*kRadTolerance*secRMin )
942       {                                           927       {
943         // At radius greater than real & imagi << 928 // At radius greater than real & imaginary cones
944         // -> 2nd root, with zi check          << 929 // -> 2nd root, with zi check
945                                                   930 
946         b = nt2/nt1 ;                             931         b = nt2/nt1 ;
947         c = nt3/nt1 ;                             932         c = nt3/nt1 ;
948         d = b*b-c ;                            << 933   d = b*b-c ;
949         if (d >= 0)   // > 0                   << 934   if (d > 0)
950         {                                      << 935   {
951            if(b>0){sd = c/( -b-std::sqrt(d));} << 936     s = -b + sqrt(d) ;
952            else   {sd = -b + std::sqrt(d) ;}   << 937 
953                                                << 938     if ( s > 0 )
954           if ( sd >= 0 )   // > 0              << 939     {
955           {                                    << 940       zi = p.z() + s*v.z() ;
956             if ( sd>dRmax ) // Avoid rounding  << 941 
957             {               // 64 bits systems << 942       if ( fabs(zi) <= tolODz )
958               G4double fTerm = sd-std::fmod(sd << 943       {
959               sd = fTerm + DistanceToIn(p+fTer << 944         if ( seg )
960             }                                  << 945         {
961             zi = p.z() + sd*v.z() ;            << 946     xi     = p.x() + s*v.x() ;
962                                                << 947     yi     = p.y() + s*v.y() ;
963             if ( std::fabs(zi) <= tolODz )     << 948     ri     = rMinAv + zi*tanRMin ;
964             {                                  << 949     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
965               if ( !fPhiFullCone )             << 950 
966               {                                << 951     if (cosPsi >= cosHDPhiOT) snxt = s ; 
967                 xi     = p.x() + sd*v.x() ;    << 952         }
968                 yi     = p.y() + sd*v.y() ;    << 953         else return s ;
969                 ri     = rMinAv + zi*tanRMin ; << 954       }
970                 cosPsi = (xi*cosCPhi + yi*sinC << 955     }
971                                                << 956   }
972                 if (cosPsi >= cosHDPhiIT)      << 
973                 {                              << 
974                   if ( sd > halfRadTolerance ) << 
975                   else                         << 
976                   {                            << 
977                     // Calculate a normal vect << 
978                                                << 
979                     risec  = std::sqrt(xi*xi + << 
980                     Normal = G4ThreeVector(-xi << 
981                     if ( Normal.dot(v) <= 0 )  << 
982                   }                            << 
983                 }                              << 
984               }                                << 
985               else                             << 
986               {                                << 
987                 if ( sd > halfRadTolerance )   << 
988                 else                           << 
989                 {                              << 
990                   // Calculate a normal vector << 
991                                                << 
992                   xi     = p.x() + sd*v.x() ;  << 
993                   yi     = p.y() + sd*v.y() ;  << 
994                   risec  = std::sqrt(xi*xi + y << 
995                   Normal = G4ThreeVector(-xi/r << 
996                   if ( Normal.dot(v) <= 0 )  { << 
997                 }                              << 
998               }                                << 
999             }                                  << 
1000           }                                   << 
1001         }                                     << 
1002       }                                          957       }
1003       else  if ( nt3 < -rin*kRadTolerance*sec    958       else  if ( nt3 < -rin*kRadTolerance*secRMin )
1004       {                                          959       {
1005         // Within radius of inner cone (real  << 960 // Within radius of inner cone (real or imaginary)
1006         // -> Try 2nd root, with checking int << 961 // -> Try 2nd root, with checking intersection is with real cone
1007         // -> If check fails, try 1st root, a << 962 // -> If check fails, try 1st root, also checking intersection is on real cone
1008         //    on real cone                    << 
1009                                                  963 
1010         b = nt2/nt1 ;                            964         b = nt2/nt1 ;
1011         c = nt3/nt1 ;                         << 965   c = nt3/nt1 ;
1012         d = b*b - c ;                         << 966   d = b*b - c ;
1013                                                  967 
1014         if ( d >= 0 )  // > 0                 << 968   if ( d > 0 )
1015         {                                        969         {
1016           if (b>0) { sd = c/(-b-std::sqrt(d)) << 970     s  = -b + sqrt(d) ;
1017           else     { sd = -b + std::sqrt(d);  << 971     zi = p.z() + s*v.z() ;
1018           zi = p.z() + sd*v.z() ;             << 972     ri = rMinAv + zi*tanRMin ;
1019           ri = rMinAv + zi*tanRMin ;          << 973 
1020                                               << 974     if ( ri >= 0 )
1021           if ( ri > 0 )                       << 975     {
1022           {                                   << 976       if ( s > 0 && fabs(zi) <= tolODz )
1023             if ( (sd >= 0) && (std::fabs(zi)  << 977       {
1024             {                                 << 978         if ( seg )
1025               if ( sd>dRmax ) // Avoid roundi << 979         {
1026               {               // seen on 64 b << 980           xi     = p.x() + s*v.x() ;
1027                 G4double fTerm = sd-std::fmod << 981     yi     = p.y() + s*v.y() ;
1028                 sd = fTerm + DistanceToIn(p+f << 982     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1029               }                               << 983 
1030               if ( !fPhiFullCone )            << 984     if (cosPsi >= cosHDPhiOT) snxt = s ; 
1031               {                               << 985         }
1032                 xi     = p.x() + sd*v.x() ;   << 986         else  return s ;
1033                 yi     = p.y() + sd*v.y() ;   << 987       }
1034                 cosPsi = (xi*cosCPhi + yi*sin << 988     }
1035                                               << 989     else
1036                 if (cosPsi >= cosHDPhiOT)     << 990     {
1037                 {                             << 991       s  = -b - sqrt(d) ;
1038                   if ( sd > halfRadTolerance  << 992       zi = p.z() + s*v.z() ;
1039                   else                        << 993       ri = rMinAv + zi*tanRMin ;
1040                   {                           << 994 
1041                     // Calculate a normal vec << 995       if ( s > 0 && ri >= 0 && fabs(zi) <= tolODz )
1042                                               << 996       {
1043                     risec  = std::sqrt(xi*xi  << 997         if ( seg )
1044                     Normal = G4ThreeVector(-x << 998         {
1045                     if ( Normal.dot(v) <= 0 ) << 999           xi     = p.x() + s*v.x() ;
1046                   }                           << 1000     yi     = p.y() + s*v.y() ;
1047                 }                             << 1001     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1048               }                               << 1002 
1049               else                            << 1003     if (cosPsi >= cosHDPhiOT) snxt = s ;
1050               {                               << 1004         }
1051                 if( sd > halfRadTolerance )   << 1005         else  return s ;
1052                 else                          << 1006       }
1053                 {                             << 1007     }
1054                   // Calculate a normal vecto << 1008   }
1055                                               << 
1056                   xi     = p.x() + sd*v.x() ; << 
1057                   yi     = p.y() + sd*v.y() ; << 
1058                   risec  = std::sqrt(xi*xi +  << 
1059                   Normal = G4ThreeVector(-xi/ << 
1060                   if ( Normal.dot(v) <= 0 )   << 
1061                 }                             << 
1062               }                               << 
1063             }                                 << 
1064           }                                   << 
1065           else                                << 
1066           {                                   << 
1067             if (b>0) { sd = -b - std::sqrt(d) << 
1068             else     { sd = c/(-b+std::sqrt(d << 
1069             zi = p.z() + sd*v.z() ;           << 
1070             ri = rMinAv + zi*tanRMin ;        << 
1071                                               << 
1072             if ( (sd >= 0) && (ri > 0) && (st << 
1073             {                                 << 
1074               if ( sd>dRmax ) // Avoid roundi << 
1075               {               // seen on 64 b << 
1076                 G4double fTerm = sd-std::fmod << 
1077                 sd = fTerm + DistanceToIn(p+f << 
1078               }                               << 
1079               if ( !fPhiFullCone )            << 
1080               {                               << 
1081                 xi     = p.x() + sd*v.x() ;   << 
1082                 yi     = p.y() + sd*v.y() ;   << 
1083                 cosPsi = (xi*cosCPhi + yi*sin << 
1084                                               << 
1085                 if (cosPsi >= cosHDPhiIT)     << 
1086                 {                             << 
1087                   if ( sd > halfRadTolerance  << 
1088                   else                        << 
1089                   {                           << 
1090                     // Calculate a normal vec << 
1091                                               << 
1092                     risec  = std::sqrt(xi*xi  << 
1093                     Normal = G4ThreeVector(-x << 
1094                     if ( Normal.dot(v) <= 0 ) << 
1095                   }                           << 
1096                 }                             << 
1097               }                               << 
1098               else                            << 
1099               {                               << 
1100                 if ( sd > halfRadTolerance )  << 
1101                 else                          << 
1102                 {                             << 
1103                   // Calculate a normal vecto << 
1104                                               << 
1105                   xi     = p.x() + sd*v.x() ; << 
1106                   yi     = p.y() + sd*v.y() ; << 
1107                   risec  = std::sqrt(xi*xi +  << 
1108                   Normal = G4ThreeVector(-xi/ << 
1109                   if ( Normal.dot(v) <= 0 )   << 
1110                 }                             << 
1111               }                               << 
1112             }                                 << 
1113           }                                   << 
1114         }                                     << 
1115       }                                          1009       }
1116       else                                       1010       else
1117       {                                          1011       {
1118         // Within kRadTol*0.5 of inner cone ( << 1012 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1119         // ----> Check not travelling through << 1013 // ----> Check not travelling through (=>0 to in)
1120         // ----> if not:                      << 1014 // ----> if not:
1121         //    -2nd root with validity check   << 1015 //    -2nd root with validity check
1122                                               << 1016 
1123         if ( std::fabs(p.z()) <= tolODz )     << 1017         if ( fabs(p.z()) <= tolODz )
1124         {                                     << 1018   {
1125           if ( nt2 > 0 )                      << 1019     if ( nt2 > 0 )
1126           {                                   << 1020     {
1127             // Inside inner real cone, headin << 1021 // Inside inner real cone, heading outwards, inside z range
1128                                               << 1022 
1129             if ( !fPhiFullCone )              << 1023       if ( seg )
1130             {                                 << 1024       {
1131               cosPsi = (p.x()*cosCPhi + p.y() << 1025         cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/sqrt(t3) ;
1132                                               << 1026 
1133               if (cosPsi >= cosHDPhiIT)  { re << 1027         if (cosPsi >= cosHDPhiIT)  return 0.0 ;
1134             }                                 << 1028       }
1135             else  { return 0.0; }             << 1029       else  return 0.0 ;
1136           }                                   << 1030     }
1137           else                                << 1031     else
1138           {                                   << 1032     {
1139             // Within z extent, but not trave << 1033 // Within z extent, but not travelling through
1140             // -> 2nd root or kInfinity if 1s << 1034 // -> 2nd root or kInfinity if 1st root on imaginary cone
1141                                               << 1035 
1142             b = nt2/nt1 ;                     << 1036       b = nt2/nt1 ;
1143             c = nt3/nt1 ;                     << 1037       c = nt3/nt1 ;
1144             d = b*b - c ;                     << 1038       d = b*b - c ;
1145                                               << 1039 
1146             if ( d >= 0 )   // > 0            << 1040       if ( d > 0 )
1147             {                                 << 1041       {
1148               if (b>0) { sd = -b - std::sqrt( << 1042         s  = -b - sqrt(d) ;
1149               else     { sd = c/(-b+std::sqrt << 1043         zi = p.z() + s*v.z() ;
1150               zi = p.z() + sd*v.z() ;         << 1044         ri = rMinAv + zi*tanRMin ;
1151               ri = rMinAv + zi*tanRMin ;      << 1045 
1152                                               << 1046         if ( ri > 0 )   // 2nd root
1153               if ( ri > 0 )   // 2nd root     << 1047         {
1154               {                               << 1048     s  = -b + sqrt(d) ;
1155                 if (b>0) { sd = c/(-b-std::sq << 1049     zi = p.z() + s*v.z() ;
1156                 else     { sd = -b + std::sqr << 1050 
1157                                               << 1051     if ( s > 0 && fabs(zi) <= tolODz )
1158                 zi = p.z() + sd*v.z() ;       << 1052     {
1159                                               << 1053       if ( seg )
1160                 if ( (sd >= 0) && (std::fabs( << 1054       {
1161                 {                             << 1055         xi     = p.x() + s*v.x() ;
1162                   if ( sd>dRmax ) // Avoid ro << 1056         yi     = p.y() + s*v.y() ;
1163                   {               // seen on  << 1057         cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1164                     G4double fTerm = sd-std:: << 1058 
1165                     sd = fTerm + DistanceToIn << 1059         if ( cosPsi >= cosHDPhiOT )  snxt = s ; 
1166                   }                           << 1060       }
1167                   if ( !fPhiFullCone )        << 1061       else return s ;
1168                   {                           << 1062     }
1169                     xi     = p.x() + sd*v.x() << 1063         }
1170                     yi     = p.y() + sd*v.y() << 1064         else  return kInfinity ;
1171                     ri     = rMinAv + zi*tanR << 1065       }
1172                     cosPsi = (xi*cosCPhi + yi << 1066     }
1173                                               << 1067   }
1174                     if ( cosPsi >= cosHDPhiIT << 
1175                   }                           << 
1176                   else  { return sd; }        << 
1177                 }                             << 
1178               }                               << 
1179               else  { return kInfinity; }     << 
1180             }                                 << 
1181           }                                   << 
1182         }                                     << 
1183         else   // 2nd root                       1068         else   // 2nd root
1184         {                                     << 1069   {
1185           b = nt2/nt1 ;                       << 1070     b = nt2/nt1 ;
1186           c = nt3/nt1 ;                       << 1071     c = nt3/nt1 ;
1187           d = b*b - c ;                       << 1072     d = b*b - c ;
1188                                               << 1073 
1189           if ( d > 0 )                        << 1074     if ( d > 0 )
1190           {                                   << 1075     { 
1191             if (b>0) { sd = c/(-b-std::sqrt(d << 1076       s  = -b + sqrt(d) ;
1192             else     { sd = -b + std::sqrt(d) << 1077       zi = p.z() + s*v.z() ;
1193             zi = p.z() + sd*v.z() ;           << 1078 
1194                                               << 1079       if ( s > 0 && fabs(zi) <= tolODz )
1195             if ( (sd >= 0) && (std::fabs(zi)  << 1080       {
1196             {                                 << 1081         if ( seg )
1197               if ( sd>dRmax ) // Avoid roundi << 1082         {
1198               {               // seen on 64 b << 1083     xi     = p.x() + s*v.x() ;
1199                 G4double fTerm = sd-std::fmod << 1084           yi     = p.y() + s*v.y() ;
1200                 sd = fTerm + DistanceToIn(p+f << 1085     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1201               }                               << 1086 
1202               if ( !fPhiFullCone )            << 1087     if (cosPsi >= cosHDPhiOT) snxt = s ; 
1203               {                               << 1088         }
1204                 xi     = p.x() + sd*v.x();    << 1089         else  return s ;
1205                 yi     = p.y() + sd*v.y();    << 1090       }
1206                 ri     = rMinAv + zi*tanRMin  << 1091     }
1207                 cosPsi = (xi*cosCPhi + yi*sin << 1092   }
1208                                               << 
1209                 if (cosPsi >= cosHDPhiIT)  {  << 
1210               }                               << 
1211               else  { return sd; }            << 
1212             }                                 << 
1213           }                                   << 
1214         }                                     << 
1215       }                                          1093       }
1216     }                                            1094     }
1217   }                                              1095   }
1218                                                  1096 
1219   // Phi segment intersection                 << 1097 // Phi segment intersection
1220   //                                          << 1098 //
1221   // o Tolerant of points inside phi planes b << 1099 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1222   //                                          << 1100 //
1223   // o NOTE: Large duplication of code betwee << 1101 // o NOTE: Large duplication of code between sphi & ephi checks
1224   //         -> only diffs: sphi -> ephi, Com << 1102 //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1225   //            intersection check <=0 -> >=0 << 1103 //            intersection check <=0 -> >=0
1226   //         -> Should use some form of loop  << 1104 //         -> Should use some form of loop Construct
1227                                                  1105 
1228   if ( !fPhiFullCone )                        << 1106   if ( seg )
1229   {                                              1107   {
1230     // First phi surface (starting phi)       << 1108 // First phi surface (`S'tarting phi)
1231                                                  1109 
                                                   >> 1110     sinSPhi = sin(fSPhi) ;
                                                   >> 1111     cosSPhi = cos(fSPhi) ;
1232     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;    1112     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
1233                                               << 1113                     
1234     if ( Comp < 0 )    // Component in outwar    1114     if ( Comp < 0 )    // Component in outwards normal dirn
1235     {                                            1115     {
1236       Dist = (p.y()*cosSPhi - p.x()*sinSPhi)     1116       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1237                                                  1117 
1238       if (Dist < halfCarTolerance)            << 1118       if (Dist < kCarTolerance*0.5)
1239       {                                          1119       {
1240         sd = Dist/Comp ;                      << 1120         s = Dist/Comp ;
1241                                               << 
1242         if ( sd < snxt )                      << 
1243         {                                     << 
1244           if ( sd < 0 )  { sd = 0.0; }        << 
1245                                               << 
1246           zi = p.z() + sd*v.z() ;             << 
1247                                               << 
1248           if ( std::fabs(zi) <= tolODz )      << 
1249           {                                   << 
1250             xi        = p.x() + sd*v.x() ;    << 
1251             yi        = p.y() + sd*v.y() ;    << 
1252             rhoi2     = xi*xi + yi*yi ;       << 
1253             tolORMin2 = (rMinOAv + zi*tanRMin << 
1254             tolORMax2 = (rMaxOAv + zi*tanRMax << 
1255                                               << 
1256             if ( (rhoi2 >= tolORMin2) && (rho << 
1257             {                                 << 
1258               // z and r intersections good - << 
1259               // correct half-plane           << 
1260                                               << 
1261               if ((yi*cosCPhi - xi*sinCPhi) < << 
1262             }                                 << 
1263           }                                   << 
1264         }                                     << 
1265       }                                       << 
1266     }                                         << 
1267                                               << 
1268     // Second phi surface (Ending phi)        << 
1269                                                  1121 
                                                   >> 1122   if ( s < snxt )
                                                   >> 1123   {
                                                   >> 1124     if ( s < 0 ) s = 0.0 ;
                                                   >> 1125 
                                                   >> 1126     zi = p.z() + s*v.z() ;
                                                   >> 1127 
                                                   >> 1128     if ( fabs(zi) <= tolODz )
                                                   >> 1129     {
                                                   >> 1130       xi        = p.x() + s*v.x() ;
                                                   >> 1131       yi        = p.y() + s*v.y() ;
                                                   >> 1132       rho2      = xi*xi + yi*yi ;
                                                   >> 1133       tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
                                                   >> 1134       tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
                                                   >> 1135 
                                                   >> 1136       if ( rho2 >= tolORMin2 && rho2 <= tolORMax2 )
                                                   >> 1137       {
                                                   >> 1138 // z and r intersections good - check intersecting with correct half-plane
                                                   >> 1139 
                                                   >> 1140         if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) snxt = s ;
                                                   >> 1141       }    
                                                   >> 1142     }
                                                   >> 1143   }
                                                   >> 1144       }
                                                   >> 1145     }    
                                                   >> 1146 // Second phi surface (`E'nding phi)
                                                   >> 1147 
                                                   >> 1148     ePhi    = fSPhi + fDPhi ;
                                                   >> 1149     sinEPhi = sin(ePhi) ;
                                                   >> 1150     cosEPhi = cos(ePhi) ;
1270     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi    1151     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1271                                               << 1152         
1272     if ( Comp < 0 )   // Component in outward    1153     if ( Comp < 0 )   // Component in outwards normal dirn
1273     {                                            1154     {
1274       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi)    1155       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1275       if (Dist < halfCarTolerance)            << 1156       if (Dist < kCarTolerance*0.5)
1276       {                                          1157       {
1277         sd = Dist/Comp ;                      << 1158         s = Dist/Comp ;
1278                                                  1159 
1279         if ( sd < snxt )                      << 1160   if ( s < snxt )
1280         {                                     << 1161   {
1281           if ( sd < 0 )  { sd = 0.0; }        << 1162     if ( s < 0 )  s = 0.0 ;
1282                                                  1163 
1283           zi = p.z() + sd*v.z() ;             << 1164     zi = p.z() + s*v.z() ;
1284                                                  1165 
1285           if (std::fabs(zi) <= tolODz)        << 1166     if (fabs(zi) <= tolODz)
1286           {                                   << 1167     {
1287             xi        = p.x() + sd*v.x() ;    << 1168       xi        = p.x() + s*v.x() ;
1288             yi        = p.y() + sd*v.y() ;    << 1169       yi        = p.y() + s*v.y() ;
1289             rhoi2     = xi*xi + yi*yi ;       << 1170       rho2      = xi*xi + yi*yi ;
1290             tolORMin2 = (rMinOAv + zi*tanRMin << 1171       tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1291             tolORMax2 = (rMaxOAv + zi*tanRMax << 1172       tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1292                                                  1173 
1293             if ( (rhoi2 >= tolORMin2) && (rho << 1174       if ( rho2 >= tolORMin2 && rho2 <= tolORMax2 )
1294             {                                 << 1175       {
1295               // z and r intersections good - << 1176 // z and r intersections good - check intersecting with correct half-plane
1296               // correct half-plane           << 
1297                                                  1177 
1298               if ( (yi*cosCPhi - xi*sinCPhi)  << 1178         if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  snxt = s ;
1299             }                                 << 1179       }    
1300           }                                   << 1180     }
1301         }                                     << 1181   }
1302       }                                          1182       }
1303     }                                            1183     }
1304   }                                              1184   }
1305   if (snxt < halfCarTolerance)  { snxt = 0.;  << 1185   if (snxt < kCarTolerance*0.5) snxt = 0.;
1306                                               << 
1307   return snxt ;                                  1186   return snxt ;
1308 }                                                1187 }
1309                                                  1188 
1310 /////////////////////////////////////////////    1189 //////////////////////////////////////////////////////////////////////////////
1311 //                                               1190 // 
1312 // Calculate distance (<= actual) to closest     1191 // Calculate distance (<= actual) to closest surface of shape from outside
1313 // - Calculate distance to z, radial planes      1192 // - Calculate distance to z, radial planes
1314 // - Only to phi planes if outside phi extent    1193 // - Only to phi planes if outside phi extent
1315 // - Return 0 if point inside                    1194 // - Return 0 if point inside
1316                                                  1195 
1317 G4double G4Cons::DistanceToIn(const G4ThreeVe    1196 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
1318 {                                                1197 {
1319   G4double safe=0.0, rho, safeR1, safeR2, saf << 1198   G4double safe, rho, safeR1, safeR2, safeZ ;
1320   G4double tanRMin, secRMin, pRMin ;             1199   G4double tanRMin, secRMin, pRMin ;
1321   G4double tanRMax, secRMax, pRMax ;             1200   G4double tanRMax, secRMax, pRMax ;
                                                   >> 1201   G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi ;
                                                   >> 1202   G4double cosPsi ;
1322                                                  1203 
1323   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y() << 1204   rho   = sqrt(p.x()*p.x() + p.y()*p.y()) ;
1324   safeZ = std::fabs(p.z()) - fDz ;            << 1205   safeZ = fabs(p.z()) - fDz ;
1325                                                  1206 
1326   if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )   << 1207   if ( fRmin1 || fRmin2 )
1327   {                                              1208   {
1328     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;        1209     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1329     secRMin = std::sqrt(1.0 + tanRMin*tanRMin << 1210     secRMin = sqrt(1.0 + tanRMin*tanRMin) ;
1330     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin    1211     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1331     safeR1  = (pRMin - rho)/secRMin ;            1212     safeR1  = (pRMin - rho)/secRMin ;
1332                                                  1213 
1333     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;        1214     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1334     secRMax = std::sqrt(1.0 + tanRMax*tanRMax << 1215     secRMax = sqrt(1.0 + tanRMax*tanRMax) ;
1335     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax    1216     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1336     safeR2  = (rho - pRMax)/secRMax ;            1217     safeR2  = (rho - pRMax)/secRMax ;
1337                                                  1218 
1338     if ( safeR1 > safeR2) { safe = safeR1; }  << 1219     if ( safeR1 > safeR2) safe = safeR1 ;
1339     else                  { safe = safeR2; }  << 1220     else                  safe = safeR2 ;
1340   }                                              1221   }
1341   else                                           1222   else
1342   {                                              1223   {
1343     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;        1224     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1344     secRMax = std::sqrt(1.0 + tanRMax*tanRMax << 1225     secRMax = sqrt(1.0 + tanRMax*tanRMax) ;
1345     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax    1226     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1346     safe    = (rho - pRMax)/secRMax ;            1227     safe    = (rho - pRMax)/secRMax ;
1347   }                                              1228   }
1348   if ( safeZ > safe )  { safe = safeZ; }      << 1229   if ( safeZ > safe ) safe = safeZ ;
1349                                                  1230 
1350   if ( !fPhiFullCone && (rho != 0.0) )        << 1231   if ( fDPhi < 2.0*M_PI && rho )
1351   {                                              1232   {
1352     // Psi=angle from central phi to point    << 1233     phiC    = fSPhi + fDPhi*0.5 ;
                                                   >> 1234     cosPhiC = cos(phiC) ;
                                                   >> 1235     sinPhiC = sin(phiC) ;
1353                                                  1236 
1354     cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/ << 1237 // Psi=angle from central phi to point
1355                                                  1238 
1356     if ( cosPsi < cosHDPhi ) // Point lies ou << 1239     cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
                                                   >> 1240 
                                                   >> 1241     if ( cosPsi < cos(fDPhi*0.5) ) // Point lies outside phi range
1357     {                                            1242     {
1358       if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= << 1243       if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0.0 )
1359       {                                          1244       {
1360         safePhi = std::fabs(p.x()*sinSPhi-p.y << 1245   safePhi = fabs(p.x()*sin(fSPhi) - p.y()*cos(fSPhi)) ;
1361       }                                          1246       }
1362       else                                       1247       else
1363       {                                          1248       {
1364         safePhi = std::fabs(p.x()*sinEPhi-p.y << 1249   ePhi    = fSPhi + fDPhi ;
                                                   >> 1250   safePhi = fabs(p.x()*sin(ePhi) - p.y()*cos(ePhi)) ;
1365       }                                          1251       }
1366       if ( safePhi > safe )  { safe = safePhi << 1252       if ( safePhi > safe ) safe = safePhi ;
1367     }                                            1253     }
1368   }                                              1254   }
1369   if ( safe < 0.0 )  { safe = 0.0; }          << 1255   if ( safe < 0.0 ) safe = 0.0 ;
1370                                                  1256 
1371   return safe ;                                  1257   return safe ;
1372 }                                                1258 }
1373                                                  1259 
1374 /////////////////////////////////////////////    1260 ///////////////////////////////////////////////////////////////
1375 //                                               1261 //
1376 // Calculate distance to surface of shape fro << 1262 // Calculate distance to surface of shape from `inside', allowing for tolerance
1377 // - Only Calc rmax intersection if no valid     1263 // - Only Calc rmax intersection if no valid rmin intersection
1378                                                  1264 
1379 G4double G4Cons::DistanceToOut( const G4Three    1265 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
1380                                 const G4Three << 1266               const G4ThreeVector& v,
1381                                 const G4bool  << 1267               const G4bool calcNorm,
1382                                       G4bool* << 1268                     G4bool *validNorm,
1383                                       G4Three << 1269                     G4ThreeVector *n) const
1384 {                                                1270 {
1385   ESide side = kNull, sider = kNull, sidephi     1271   ESide side = kNull, sider = kNull, sidephi = kNull;
1386                                                  1272 
1387   G4double snxt,srd,sphi,pdist ;              << 1273   G4double snxt,sr,sphi,pdist ;
1388                                                  1274 
1389   G4double tanRMax, secRMax, rMaxAv ;  // Dat << 1275   G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1390   G4double tanRMin, secRMin, rMinAv ;  // Dat << 1276   G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1391                                                  1277 
1392   G4double t1, t2, t3, rout, rin, nt1, nt2, n    1278   G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1393   G4double b, c, d, sr2, sr3 ;                   1279   G4double b, c, d, sr2, sr3 ;
1394                                                  1280 
1395   // Vars for intersection within tolerance   << 1281 // Vars for intersection within tolerance
1396                                                  1282 
1397   ESide sidetol = kNull ;                     << 1283   ESide    sidetol ;
1398   G4double slentol = kInfinity ;                 1284   G4double slentol = kInfinity ;
1399                                                  1285 
1400   // Vars for phi intersection:               << 1286 // Vars for phi intersection:
1401                                                  1287 
                                                   >> 1288   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
                                                   >> 1289   G4double cPhi, sinCPhi, cosCPhi ;
1402   G4double pDistS, compS, pDistE, compE, sphi    1290   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1403   G4double zi, ri, deltaRoi2 ;                   1291   G4double zi, ri, deltaRoi2 ;
1404                                                  1292 
1405   // Z plane intersection                     << 1293 // Z plane intersection
1406                                                  1294 
1407   if ( v.z() > 0.0 )                             1295   if ( v.z() > 0.0 )
1408   {                                              1296   {
1409     pdist = fDz - p.z() ;                        1297     pdist = fDz - p.z() ;
1410                                                  1298 
1411     if (pdist > halfCarTolerance)             << 1299     if (pdist > kCarTolerance*0.5)
1412     {                                            1300     {
1413       snxt = pdist/v.z() ;                       1301       snxt = pdist/v.z() ;
1414       side = kPZ ;                               1302       side = kPZ ;
1415     }                                            1303     }
1416     else                                         1304     else
1417     {                                            1305     {
1418       if (calcNorm)                              1306       if (calcNorm)
1419       {                                          1307       {
1420         *n         = G4ThreeVector(0,0,1) ;   << 1308   *n         = G4ThreeVector(0,0,1) ;
1421         *validNorm = true ;                   << 1309   *validNorm = true ;
1422       }                                          1310       }
1423       return  snxt = 0.0;                     << 1311       return snxt = 0.0 ;
1424     }                                            1312     }
1425   }                                              1313   }
1426   else if ( v.z() < 0.0 )                        1314   else if ( v.z() < 0.0 )
1427   {                                              1315   {
1428     pdist = fDz + p.z() ;                        1316     pdist = fDz + p.z() ;
1429                                                  1317 
1430     if ( pdist > halfCarTolerance)            << 1318     if ( pdist > kCarTolerance*0.5)
1431     {                                            1319     {
1432       snxt = -pdist/v.z() ;                      1320       snxt = -pdist/v.z() ;
1433       side = kMZ ;                               1321       side = kMZ ;
1434     }                                            1322     }
1435     else                                         1323     else
1436     {                                            1324     {
1437       if ( calcNorm )                            1325       if ( calcNorm )
1438       {                                          1326       {
1439         *n         = G4ThreeVector(0,0,-1) ;  << 1327   *n         = G4ThreeVector(0,0,-1) ;
1440         *validNorm = true ;                   << 1328   *validNorm = true ;
1441       }                                          1329       }
1442       return snxt = 0.0 ;                        1330       return snxt = 0.0 ;
1443     }                                            1331     }
1444   }                                              1332   }
1445   else     // Travel perpendicular to z axis     1333   else     // Travel perpendicular to z axis
1446   {                                              1334   {
1447     snxt = kInfinity ;                        << 1335     snxt = kInfinity ;    
1448     side = kNull ;                               1336     side = kNull ;
1449   }                                              1337   }
1450                                                  1338 
1451   // Radial Intersections                     << 1339 // Radial Intersections
1452   //                                          << 1340 //
1453   // Intersection with outer cone (possible r << 1341 // Intersection with outer cone (possible return) and
1454   //                   inner cone (must also  << 1342 //                   inner cone (must also check phi)
1455   //                                          << 1343 //
1456   // Intersection point (xi,yi,zi) on line x= << 1344 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1457   //                                          << 1345 //
1458   // Intersects with x^2+y^2=(a*z+b)^2        << 1346 // Intersects with x^2+y^2=(a*z+b)^2
1459   //                                          << 1347 //
1460   // where a=tanRMax or tanRMin               << 1348 // where a=tanRMax or tanRMin
1461   //       b=rMaxAv  or rMinAv                << 1349 //       b=rMaxAv  or rMinAv
1462   //                                          << 1350 //
1463   // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*v << 1351 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
1464   //     t1                        t2         << 1352 //     t1                        t2                      t3  
1465   //                                          << 1353 //
1466   //  \--------u-------/       \-----------v- << 1354 //  \--------u-------/       \-----------v----------/ \---------w--------/
1467                                                  1355 
1468   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;          1356   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1469   secRMax = std::sqrt(1.0 + tanRMax*tanRMax)  << 1357   secRMax = sqrt(1.0 + tanRMax*tanRMax) ;
1470   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;              1358   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
1471                                                  1359 
1472                                                  1360 
1473   t1   = 1.0 - v.z()*v.z() ;      // since v  << 1361   t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1474   t2   = p.x()*v.x() + p.y()*v.y() ;             1362   t2   = p.x()*v.x() + p.y()*v.y() ;
1475   t3   = p.x()*p.x() + p.y()*p.y() ;             1363   t3   = p.x()*p.x() + p.y()*p.y() ;
1476   rout = tanRMax*p.z() + rMaxAv ;                1364   rout = tanRMax*p.z() + rMaxAv ;
1477                                                  1365 
1478   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z())     1366   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1479   nt2 = t2 - tanRMax*v.z()*rout ;                1367   nt2 = t2 - tanRMax*v.z()*rout ;
1480   nt3 = t3 - rout*rout ;                         1368   nt3 = t3 - rout*rout ;
1481                                                  1369 
1482   if (v.z() > 0.0)                               1370   if (v.z() > 0.0)
1483   {                                              1371   {
1484     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 << 1372     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1485                 - fRmax2*(fRmax2 + kRadTolera << 
1486   }                                              1373   }
1487   else if (v.z() < 0.0)                       << 1374   else if ( v.z() < 0.0 )
1488   {                                           << 
1489     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 << 
1490                 - fRmax1*(fRmax1 + kRadTolera << 
1491   }                                           << 
1492   else                                        << 
1493   {                                              1375   {
1494     deltaRoi2 = 1.0;                          << 1376     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1495   }                                              1377   }
                                                   >> 1378   else deltaRoi2 = 1.0 ;
1496                                                  1379 
1497   if ( (nt1 != 0.0) && (deltaRoi2 > 0.0) )    << 1380   if ( nt1 && deltaRoi2 > 0.0 )  
1498   {                                              1381   {
1499     // Equation quadratic => 2 roots : second << 1382 // Equation quadratic => 2 roots : second root must be leaving
1500                                                  1383 
1501     b = nt2/nt1 ;                                1384     b = nt2/nt1 ;
1502     c = nt3/nt1 ;                                1385     c = nt3/nt1 ;
1503     d = b*b - c ;                                1386     d = b*b - c ;
1504                                                  1387 
1505     if ( d >= 0 )                                1388     if ( d >= 0 )
1506     {                                            1389     {
1507       // Check if on outer cone & heading out << 1390 // Check if on outer cone & heading outwards
1508       // NOTE: Should use rho-rout>-kRadToler << 1391 // NOTE: Should use rho-rout>-kRadtolerance*0.5
1509                                               << 1392         
1510       if (nt3 > -halfRadTolerance && nt2 >= 0 << 1393       if (nt3 > -kRadTolerance*0.5 && nt2 >= 0 )
1511       {                                       << 1394       {
1512         if (calcNorm)                         << 1395   if (calcNorm)
1513         {                                     << 1396   {
1514           risec      = std::sqrt(t3)*secRMax  << 1397     risec      = sqrt(t3)*secRMax ;
1515           *validNorm = true ;                 << 1398     *validNorm = true ;
1516           *n         = G4ThreeVector(p.x()/ri << 1399     *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
1517         }                                     << 1400   }
1518         return snxt=0 ;                       << 1401   return snxt=0 ;
1519       }                                          1402       }
1520       else                                       1403       else
1521       {                                          1404       {
1522         sider = kRMax  ;                      << 1405   sider = kRMax  ;
1523         if (b>0) { srd = -b - std::sqrt(d);   << 1406   sr    = -b - sqrt(d) ; // was +srqrt(d), vmg 28.04.99
1524         else     { srd = c/(-b+std::sqrt(d))  << 1407   zi    = p.z() + sr*v.z() ;
1525                                               << 1408   ri    = tanRMax*zi + rMaxAv ;
1526         zi    = p.z() + srd*v.z() ;           << 1409           
1527         ri    = tanRMax*zi + rMaxAv ;         << 1410   if ( (ri >= 0) && (-kRadTolerance*0.5 <= sr) && 
1528                                               << 1411                           ( sr <= kRadTolerance*0.5)     )
1529         if ((ri >= 0) && (-halfRadTolerance < << 1412   {
1530         {                                     << 1413 // An intersection within the tolerance
1531           // An intersection within the toler << 1414 //   we will Store it in case it is good -
1532           //   we will Store it in case it is << 1415 // 
1533           //                                  << 1416     slentol = sr ;
1534           slentol = srd ;                     << 1417     sidetol = kRMax ;
1535           sidetol = kRMax ;                   << 1418   }           
1536         }                                     << 1419   if ( (ri < 0) || (sr < kRadTolerance*0.5) )
1537         if ( (ri < 0) || (srd < halfRadTolera << 1420   {
1538         {                                     << 1421 // Safety: if both roots -ve ensure that sr cannot `win' distance to out
1539           // Safety: if both roots -ve ensure << 1422 
1540           //         distance to out          << 1423     sr2 = -b + sqrt(d) ;
1541                                               << 1424     zi  = p.z() + sr2*v.z() ;
1542           if (b>0) { sr2 = c/(-b-std::sqrt(d) << 1425     ri  = tanRMax*zi + rMaxAv ;
1543           else     { sr2 = -b + std::sqrt(d); << 1426 
1544           zi  = p.z() + sr2*v.z() ;           << 1427     if (ri >= 0 && sr2 > kRadTolerance*0.5)  sr = sr2 ;
1545           ri  = tanRMax*zi + rMaxAv ;         << 1428     else
1546                                               << 1429     {
1547           if ((ri >= 0) && (sr2 > halfRadTole << 1430       sr = kInfinity ;
1548           {                                   << 1431 
1549             srd = sr2;                        << 1432       if(    (-kRadTolerance*0.5 <= sr2)
1550           }                                   << 1433            && ( sr2 <= kRadTolerance*0.5)  )
1551           else                                << 1434       {
1552           {                                   << 1435 // An intersection within the tolerance. Storing it in case it is good.
1553             srd = kInfinity ;                 << 1436 
1554                                               << 1437         slentol = sr2 ;
1555             if( (-halfRadTolerance <= sr2) && << 1438         sidetol = kRMax ;
1556             {                                 << 1439       }
1557               // An intersection within the t << 1440     }
1558               // Storing it in case it is goo << 1441   }
1559                                               << 
1560               slentol = sr2 ;                 << 
1561               sidetol = kRMax ;               << 
1562             }                                 << 
1563           }                                   << 
1564         }                                     << 
1565       }                                          1442       }
1566     }                                            1443     }
1567     else                                         1444     else
1568     {                                            1445     {
1569       // No intersection with outer cone & no << 1446 // No intersection with outer cone & not parallel -> already outside, no
1570       // -> already outside, no intersection  << 1447 // intersection
1571                                                  1448 
1572       if ( calcNorm )                            1449       if ( calcNorm )
1573       {                                          1450       {
1574         risec      = std::sqrt(t3)*secRMax;   << 1451   risec      = sqrt(t3)*secRMax ;
1575         *validNorm = true;                    << 1452   *validNorm = true ;
1576         *n         = G4ThreeVector(p.x()/rise << 1453   *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
1577       }                                          1454       }
1578       return snxt = 0.0 ;                        1455       return snxt = 0.0 ;
1579     }                                            1456     }
1580   }                                              1457   }
1581   else if ( (nt2 != 0.0) && (deltaRoi2 > 0.0) << 1458   else if ( nt2 && deltaRoi2 > 0.0 )
1582   {                                              1459   {
1583     // Linear case (only one intersection) => << 1460 // Linear case (only one intersection) => point outside outer cone
1584                                                  1461 
1585     if ( calcNorm )                              1462     if ( calcNorm )
1586     {                                            1463     {
1587       risec      = std::sqrt(t3)*secRMax;     << 1464       risec      = sqrt(t3)*secRMax ;
1588       *validNorm = true;                      << 1465       *validNorm = true ;
1589       *n         = G4ThreeVector(p.x()/risec, << 1466       *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
1590     }                                            1467     }
1591     return snxt = 0.0 ;                          1468     return snxt = 0.0 ;
1592   }                                              1469   }
1593   else                                           1470   else
1594   {                                              1471   {
1595     // No intersection -> parallel to outer c << 1472 // No intersection -> parallel to outer cone => Z or inner cone intersection
1596     // => Z or inner cone intersection        << 
1597                                                  1473 
1598     srd = kInfinity ;                         << 1474        sr = kInfinity ;
1599   }                                              1475   }
1600                                                  1476 
1601   // Check possible intersection within toler << 1477 // Check possible intersection within tolerance
1602                                                  1478 
1603   if ( slentol <= halfCarTolerance )          << 1479   if ( slentol <= kCarTolerance*0.5 )
1604   {                                              1480   {
1605     // An intersection within the tolerance w << 1481 // An intersection within the tolerance was found.  
1606     // We must accept it only if the momentum << 1482 // We must accept it only if the momentum points outwards.  
1607     //                                        << 1483 //
1608     // G4ThreeVector ptTol ;  // The point of << 1484 // G4ThreeVector ptTol ;  // The point of the intersection  
1609     // ptTol= p + slentol*v ;                 << 1485 // ptTol= p + slentol*v ;
1610     // ri=tanRMax*zi+rMaxAv ;                 << 1486 // ri=tanRMax*zi+rMaxAv ;
1611     //                                        << 1487 //
1612     // Calculate a normal vector,  as below   << 1488 // Calculate a normal vector,  as below
1613                                                  1489 
1614     xi    = p.x() + slentol*v.x();            << 1490     xi    = p.x() + slentol*v.x() ;
1615     yi    = p.y() + slentol*v.y();            << 1491     yi    = p.y() + slentol*v.y() ;
1616     risec = std::sqrt(xi*xi + yi*yi)*secRMax; << 1492     risec = sqrt(xi*xi + yi*yi)*secRMax ;
1617     G4ThreeVector Normal = G4ThreeVector(xi/r << 1493     G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1618                                                  1494 
1619     if ( Normal.dot(v) > 0 )    // We will le    1495     if ( Normal.dot(v) > 0 )    // We will leave the Cone immediatelly
1620     {                                            1496     {
1621       if ( calcNorm )                            1497       if ( calcNorm ) 
1622       {                                          1498       {
1623         *n         = Normal.unit() ;          << 1499   *n         = Normal.unit() ;
1624         *validNorm = true ;                   << 1500   *validNorm = true ;
1625       }                                          1501       }
1626       return snxt = 0.0 ;                        1502       return snxt = 0.0 ;
1627     }                                            1503     }
1628     else // On the surface, but not heading o    1504     else // On the surface, but not heading out so we ignore this intersection
1629     {    //                                   << 1505     {    //                                    (as it is within tolerance).  
1630       slentol = kInfinity ;                      1506       slentol = kInfinity ;
1631     }                                            1507     }
1632   }                                              1508   }
1633                                                  1509 
1634   // Inner Cone intersection                  << 1510 // Inner Cone intersection
1635                                                  1511 
1636   if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )   << 1512   if ( fRmin1 || fRmin2 )
1637   {                                              1513   {
1638     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;        1514     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1639     nt1     = t1 - (tanRMin*v.z())*(tanRMin*v    1515     nt1     = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1640                                                  1516 
1641     if ( nt1 != 0.0 )                         << 1517     if ( nt1 )
1642     {                                            1518     {
1643       secRMin = std::sqrt(1.0 + tanRMin*tanRM << 1519       secRMin = sqrt(1.0 + tanRMin*tanRMin) ;
1644       rMinAv  = (fRmin1 + fRmin2)*0.5 ;          1520       rMinAv  = (fRmin1 + fRmin2)*0.5 ;    
1645       rin     = tanRMin*p.z() + rMinAv ;         1521       rin     = tanRMin*p.z() + rMinAv ;
1646       nt2     = t2 - tanRMin*v.z()*rin ;         1522       nt2     = t2 - tanRMin*v.z()*rin ;
1647       nt3     = t3 - rin*rin ;                   1523       nt3     = t3 - rin*rin ;
1648                                               << 1524       
1649       // Equation quadratic => 2 roots : firs << 1525 // Equation quadratic => 2 roots : first root must be leaving
1650                                                  1526 
1651       b = nt2/nt1 ;                              1527       b = nt2/nt1 ;
1652       c = nt3/nt1 ;                              1528       c = nt3/nt1 ;
1653       d = b*b - c ;                              1529       d = b*b - c ;
1654                                                  1530 
1655       if ( d >= 0.0 )                         << 1531       if (d >= 0.0 )
1656       {                                          1532       {
1657         // NOTE: should be rho-rin<kRadTolera << 1533 // NOTE: should be rho-rin<kRadTolerance*0.5, but using squared versions
1658         //       but using squared versions f << 1534 // for efficiency
1659                                                  1535 
1660         if (nt3 < kRadTolerance*(rin + kRadTo    1536         if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25)) 
1661         {                                     << 1537   {
1662           if ( nt2 < 0.0 )                    << 1538     if ( nt2 < 0.0 )
1663           {                                   << 1539     {
1664             if (calcNorm)  { *validNorm = fal << 1540       if (calcNorm)  *validNorm = false ;
1665             return          snxt      = 0.0;  << 1541       return          snxt      = 0.0 ;
1666           }                                   << 1542     }
1667         }                                     << 1543   }
1668         else                                  << 1544   else
1669         {                                     << 1545   {
1670           if (b>0) { sr2 = -b - std::sqrt(d); << 1546     sr2 = -b - sqrt(d) ;
1671           else     { sr2 = c/(-b+std::sqrt(d) << 1547     zi  = p.z() + sr2*v.z() ;
1672           zi  = p.z() + sr2*v.z() ;           << 1548     ri  = tanRMin*zi + rMinAv ;
1673           ri  = tanRMin*zi + rMinAv ;         << 1549 
1674                                               << 1550     if( (ri >= 0.0) && (-kRadTolerance*0.5 <= sr2) && 
1675           if( (ri>=0.0)&&(-halfRadTolerance<= << 1551                              ( sr2 <= kRadTolerance*0.5)      )
1676           {                                   << 1552     {
1677             // An intersection within the tol << 1553 // An intersection within the tolerance  storing it in case it is good.
1678             // storing it in case it is good. << 1554 
1679                                               << 1555       slentol = sr2 ;
1680             slentol = sr2 ;                   << 1556       sidetol = kRMax ;
1681             sidetol = kRMax ;                 << 1557     }
1682           }                                   << 1558     if( (ri<0) || (sr2 < kRadTolerance*0.5) )
1683           if( (ri<0) || (sr2 < halfRadToleran << 1559     {
1684           {                                   << 1560       sr3 = -b + sqrt(d) ;
1685             if (b>0) { sr3 = c/(-b-std::sqrt( << 1561 
1686             else     { sr3 = -b + std::sqrt(d << 1562 // Safety: if both roots -ve ensure that sr cannot `win' distancetoout
1687                                               << 1563 
1688             // Safety: if both roots -ve ensu << 1564       if  ( sr3 > kCarTolerance*0.5 )
1689             //         distancetoout          << 1565       {
1690                                               << 1566         if( sr3 < sr )
1691             if  ( sr3 > halfRadTolerance )    << 1567         {
1692             {                                 << 1568           zi = p.z() + sr3*v.z() ;
1693               if( sr3 < srd )                 << 1569     ri = tanRMin*zi + rMinAv ;
1694               {                               << 1570 
1695                 zi = p.z() + sr3*v.z() ;      << 1571     if ( ri >= 0.0 )
1696                 ri = tanRMin*zi + rMinAv ;    << 1572     {
1697                                               << 1573       sr=sr3 ;
1698                 if ( ri >= 0.0 )              << 1574       sider=kRMin ;
1699                 {                             << 1575           }
1700                   srd=sr3 ;                   << 1576         } 
1701                   sider=kRMin ;               << 1577       }
1702                 }                             << 1578       else if ( sr3 > -kCarTolerance*0.5 )
1703               }                               << 1579       {
1704             }                                 << 1580 // Intersection in tolerance. Store to check if it's good
1705             else if ( sr3 > -halfRadTolerance << 1581 
1706             {                                 << 1582         slentol = sr3 ;
1707               // Intersection in tolerance. S << 1583         sidetol = kRMin ;
1708                                               << 1584       }
1709               slentol = sr3 ;                 << 1585     }
1710               sidetol = kRMin ;               << 1586     else if ( sr2 < sr && sr2 > kCarTolerance*0.5 )
1711             }                                 << 1587     {
1712           }                                   << 1588       sr    = sr2 ;
1713           else if ( (sr2 < srd) && (sr2 > hal << 1589       sider = kRMin ;
1714           {                                   << 1590     }
1715             srd   = sr2 ;                     << 1591     else if (sr2 > -kCarTolerance*0.5)
1716             sider = kRMin ;                   << 1592     {
1717           }                                   << 1593 // Intersection in tolerance. Store to check if it's good
1718           else if (sr2 > -halfCarTolerance)   << 1594 
1719           {                                   << 1595       slentol = sr2 ;
1720             // Intersection in tolerance. Sto << 1596       sidetol = kRMin ;
1721                                               << 1597     }   
1722             slentol = sr2 ;                   << 1598     if( slentol <= kCarTolerance*0.5  )
1723             sidetol = kRMin ;                 << 1599     {
1724           }                                   << 1600 // An intersection within the tolerance was found. 
1725           if( slentol <= halfCarTolerance  )  << 1601 // We must accept it only if  the momentum points outwards. 
1726           {                                   << 1602 
1727             // An intersection within the tol << 1603       G4ThreeVector Normal ; 
1728             // We must accept it only if  the << 1604             
1729                                               << 1605 // Calculate a normal vector,  as below
1730             G4ThreeVector Normal ;            << 1606 
1731                                               << 1607       xi     = p.x() + slentol*v.x() ;
1732             // Calculate a normal vector,  as << 1608       yi     = p.y() + slentol*v.y() ;
1733                                               << 1609       risec  = sqrt(xi*xi + yi*yi)*secRMin ;
1734             xi     = p.x() + slentol*v.x() ;  << 1610       Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMin/secRMin) ;
1735             yi     = p.y() + slentol*v.y() ;  << 1611             
1736             if( sidetol==kRMax )              << 1612       if( Normal.dot(v) > 0 )
1737             {                                 << 1613       {
1738               risec  = std::sqrt(xi*xi + yi*y << 1614 // We will leave the Cone immediatelly
1739               Normal = G4ThreeVector(xi/risec << 1615         if( calcNorm ) 
1740             }                                 << 1616         {
1741             else                              << 1617     *n         = Normal.unit() ;
1742             {                                 << 1618     *validNorm = true ;
1743               risec  = std::sqrt(xi*xi + yi*y << 1619         }
1744               Normal = G4ThreeVector(-xi/rise << 1620         return snxt = 0.0 ;
1745             }                                 << 1621       }
1746             if( Normal.dot(v) > 0 )           << 1622       else 
1747             {                                 << 1623       { 
1748               // We will leave the cone immed << 1624 // On the surface, but not heading out so we ignore this
1749                                               << 1625 // intersection (as it is within tolerance). 
1750               if( calcNorm )                  << 1626 
1751               {                               << 1627         slentol = kInfinity ;
1752                 *n         = Normal.unit() ;  << 1628       }       
1753                 *validNorm = true ;           << 1629     }
1754               }                               << 1630   }
1755               return snxt = 0.0 ;             << 1631       }
1756             }                                 << 1632     }
1757             else                              << 1633   }
1758             {                                 << 1634 
1759               // On the surface, but not head << 1635 // Linear case => point outside inner cone ---> outer cone intersect
1760               // intersection (as it is withi << 1636 //
1761                                               << 1637 // Phi Intersection
1762               slentol = kInfinity ;           << 
1763             }                                 << 
1764           }                                   << 
1765         }                                     << 
1766       }                                       << 
1767     }                                         << 
1768   }                                           << 
1769                                               << 
1770   // Linear case => point outside inner cone  << 
1771   //                                          << 
1772   // Phi Intersection                         << 
1773                                                  1638   
1774   if ( !fPhiFullCone )                        << 1639   if ( fDPhi < 2.0*M_PI )
1775   {                                              1640   {
1776     // add angle calculation with correction  << 1641     sinSPhi = sin(fSPhi) ;
1777     // of the difference in domain of atan2 a << 1642     cosSPhi = cos(fSPhi) ;
1778                                               << 1643     ePhi    = fSPhi + fDPhi ;
1779     vphi = std::atan2(v.y(),v.x()) ;          << 1644     sinEPhi = sin(ePhi) ;
1780                                               << 1645     cosEPhi = cos(ePhi) ;
1781     if ( vphi < fSPhi - halfAngTolerance  )   << 1646     cPhi    = fSPhi + fDPhi*0.5 ;
1782     else if ( vphi > fSPhi + fDPhi + halfAngT << 1647     sinCPhi = sin(cPhi) ;
                                                   >> 1648     cosCPhi = cos(cPhi) ;
1783                                                  1649 
1784     if ( (p.x() != 0.0) || (p.y() != 0.0) )   << 1650     if ( p.x() || p.y() )   // Check if on z axis (rho not needed later)
1785     {                                            1651     {
1786       // pDist -ve when inside                << 1652 // pDist -ve when inside
1787                                                  1653 
1788       pDistS = p.x()*sinSPhi - p.y()*cosSPhi     1654       pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1789       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi    1655       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1790                                                  1656 
1791       // Comp -ve when in direction of outwar << 1657 // Comp -ve when in direction of outwards normal
1792                                                  1658 
1793       compS = -sinSPhi*v.x() + cosSPhi*v.y()     1659       compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1794       compE = sinEPhi*v.x() - cosEPhi*v.y() ;    1660       compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1795                                                  1661 
1796       sidephi = kNull ;                          1662       sidephi = kNull ;
1797                                                  1663 
1798       if( ( (fDPhi <= pi) && ( (pDistS <= hal << 1664       if (pDistS <= 0 && pDistE <= 0)
1799                             && (pDistE <= hal << 
1800          || ( (fDPhi >  pi) && ((pDistS <=  h << 
1801                               || (pDistE <=   << 
1802       {                                          1665       {
1803         // Inside both phi *full* planes      << 1666 // Inside both phi *full* planes
1804         if ( compS < 0 )                      << 
1805         {                                     << 
1806           sphi = pDistS/compS ;               << 
1807           if (sphi >= -halfCarTolerance)      << 
1808           {                                   << 
1809             xi = p.x() + sphi*v.x() ;         << 
1810             yi = p.y() + sphi*v.y() ;         << 
1811                                               << 
1812             // Check intersecting with correc << 
1813             // (if not -> no intersect)       << 
1814             //                                << 
1815             if ( (std::fabs(xi)<=kCarToleranc << 
1816               && (std::fabs(yi)<=kCarToleranc << 
1817             {                                 << 
1818               sidephi= kSPhi;                 << 
1819               if ( ( fSPhi-halfAngTolerance < << 
1820                 && ( fSPhi+fDPhi+halfAngToler << 
1821               {                               << 
1822                 sphi = kInfinity;             << 
1823               }                               << 
1824             }                                 << 
1825             else                              << 
1826             if ( (yi*cosCPhi-xi*sinCPhi)>=0 ) << 
1827             {                                 << 
1828               sphi = kInfinity ;              << 
1829             }                                 << 
1830             else                              << 
1831             {                                 << 
1832               sidephi = kSPhi ;               << 
1833               if ( pDistS > -halfCarTolerance << 
1834               {                               << 
1835                 sphi = 0.0 ; // Leave by sphi << 
1836               }                               << 
1837             }                                 << 
1838           }                                   << 
1839           else                                << 
1840           {                                   << 
1841             sphi = kInfinity ;                << 
1842           }                                   << 
1843         }                                     << 
1844         else                                  << 
1845         {                                     << 
1846           sphi = kInfinity ;                  << 
1847         }                                     << 
1848                                                  1667 
1849         if ( compE < 0 )                      << 1668         if (compS < 0)
1850         {                                        1669         {
1851           sphi2 = pDistE/compE ;              << 1670     sphi = pDistS/compS ;
1852                                               << 1671     xi   = p.x() + sphi*v.x() ;
1853           // Only check further if < starting << 1672     yi   = p.y() + sphi*v.y() ;
1854           //                                  << 1673 
1855           if ( (sphi2 > -halfCarTolerance) && << 1674 // Check intersecting with correct half-plane (if not -> no intersect)
1856           {                                   << 1675 
1857             xi = p.x() + sphi2*v.x() ;        << 1676     if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) sphi = kInfinity ;
1858             yi = p.y() + sphi2*v.y() ;        << 1677     else
1859                                               << 1678     {
1860             // Check intersecting with correc << 1679       sidephi = kSPhi ;
1861                                               << 1680 
1862             if ( (std::fabs(xi)<=kCarToleranc << 1681       if ( pDistS > -kCarTolerance*0.5 ) // Leave by sphi immediately
1863               && (std::fabs(yi)<=kCarToleranc << 1682       {
1864             {                                 << 1683         sphi = 0.0 ;
1865               // Leaving via ending phi       << 1684       }     
1866                                               << 1685     }
1867               if( (fSPhi-halfAngTolerance > v << 1686   }
1868                  || (fSPhi+fDPhi+halfAngToler << 1687   else sphi = kInfinity ;
1869               {                               << 1688           
1870                 sidephi = kEPhi ;             << 1689   if ( compE < 0.0 )
1871                 if ( pDistE <= -halfCarTolera << 1690   {
1872                 else                          << 1691 // Only check further if < starting phi intersection
1873               }                               << 1692 
1874             }                                 << 1693     sphi2 = pDistE/compE ;
1875             else // Check intersecting with c << 1694 
1876             if ( yi*cosCPhi-xi*sinCPhi >= 0 ) << 1695     if ( sphi2 < sphi )
1877             {                                 << 1696     {
1878               // Leaving via ending phi       << 1697       xi = p.x() + sphi2*v.x() ;
1879                                               << 1698       yi = p.y() + sphi2*v.y() ;
1880               sidephi = kEPhi ;               << 1699 
1881               if ( pDistE <= -halfCarToleranc << 1700 // Check intersecting with correct half-plane 
1882               else                            << 1701 
1883             }                                 << 1702       if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  // Leaving via ending phi
1884           }                                   << 1703       {
1885         }                                     << 1704         sidephi = kEPhi ;
                                                   >> 1705 
                                                   >> 1706         if ( pDistE <= -kCarTolerance*0.5 )  sphi = sphi2 ;
                                                   >> 1707         else                                 sphi = 0.0   ;
                                                   >> 1708       }
                                                   >> 1709     }
                                                   >> 1710   }       
                                                   >> 1711       }
                                                   >> 1712       else if ( pDistS >= 0 && pDistE >= 0 )
                                                   >> 1713       {
                                                   >> 1714 // Outside both *full* phi planes
                                                   >> 1715 
                                                   >> 1716         if ( pDistS <= pDistE )  sidephi = kSPhi ;
                                                   >> 1717         else                     sidephi = kEPhi ;
                                                   >> 1718 
                                                   >> 1719         if ( fDPhi > M_PI )
                                                   >> 1720   {
                                                   >> 1721     if (compS < 0 && compE < 0) sphi = 0.0 ;
                                                   >> 1722     else                        sphi = kInfinity ;
                                                   >> 1723   }
                                                   >> 1724   else
                                                   >> 1725   {
                                                   >> 1726 // if towards both >=0 then once inside (after error) will remain inside
                                                   >> 1727 
                                                   >> 1728           if ( compS >= 0.0 && compE >= 0.0 ) sphi = kInfinity ;
                                                   >> 1729     else                                sphi = 0.0 ;
                                                   >> 1730   }       
                                                   >> 1731       }
                                                   >> 1732       else if ( pDistS > 0.0 && pDistE < 0.0 )
                                                   >> 1733       {
                                                   >> 1734 // Outside full starting plane, inside full ending plane
                                                   >> 1735 
                                                   >> 1736         if ( fDPhi > M_PI )
                                                   >> 1737         {
                                                   >> 1738     if ( compE < 0.0 )
                                                   >> 1739     {
                                                   >> 1740       sphi = pDistE/compE ;
                                                   >> 1741       xi   = p.x() + sphi*v.x() ;
                                                   >> 1742       yi   = p.y() + sphi*v.y() ;
                                                   >> 1743 
                                                   >> 1744 // Check intersection in correct half-plane (if not -> not leaving phi extent)
                                                   >> 1745 
                                                   >> 1746       if ( (yi*cosCPhi - xi*sinCPhi) <= 0.0 )
                                                   >> 1747       {
                                                   >> 1748         sphi = kInfinity ;
                                                   >> 1749       }
                                                   >> 1750       else    // Leaving via Ending phi
                                                   >> 1751       {
                                                   >> 1752               sidephi = kEPhi  ;
                                                   >> 1753 
                                                   >> 1754         if ( pDistE > -kCarTolerance*0.5 ) sphi = 0.0 ;
                                                   >> 1755       }
                                                   >> 1756     }
                                                   >> 1757     else sphi = kInfinity ;
                                                   >> 1758   }
                                                   >> 1759   else  // fDPhi <= pi
                                                   >> 1760   {
                                                   >> 1761     if (compS >= 0.0 )
                                                   >> 1762     {
                                                   >> 1763       if ( compE < 0.0 )
                                                   >> 1764       {           
                                                   >> 1765         sphi = pDistE/compE ;
                                                   >> 1766         xi   = p.x() + sphi*v.x() ;
                                                   >> 1767         yi   = p.y() + sphi*v.y() ;
                                                   >> 1768 
                                                   >> 1769 // Check intersection in correct half-plane (if not -> remain in extent)
                                                   >> 1770 
                                                   >> 1771         if ( (yi*cosCPhi - xi*sinCPhi) <= 0.0 )
                                                   >> 1772         {
                                                   >> 1773     sphi = kInfinity ;
                                                   >> 1774         }
                                                   >> 1775         else  // otherwise leaving via Ending phi
                                                   >> 1776         {
                                                   >> 1777     sidephi = kEPhi ;
                                                   >> 1778         }
                                                   >> 1779       }
                                                   >> 1780       else sphi = kInfinity ;
                                                   >> 1781     }
                                                   >> 1782     else  // leaving immediately by starting phi
                                                   >> 1783     {
                                                   >> 1784       sidephi = kSPhi ;
                                                   >> 1785       sphi    = 0.0 ;
                                                   >> 1786     }
                                                   >> 1787   }
1886       }                                          1788       }
1887       else                                       1789       else
1888       {                                          1790       {
1889         sphi = kInfinity ;                    << 1791 // Must be pDistS<0&&pDistE>0
1890       }                                       << 1792 // Inside full starting plane, outside full ending plane
                                                   >> 1793 
                                                   >> 1794         if (fDPhi > M_PI)
                                                   >> 1795   {
                                                   >> 1796     if ( compS < 0.0 )
                                                   >> 1797     {
                                                   >> 1798       sphi = pDistS/compS ;
                                                   >> 1799       xi   = p.x() + sphi*v.x() ;
                                                   >> 1800       yi   = p.y() + sphi*v.y() ;
                                                   >> 1801 
                                                   >> 1802 // Check intersection in correct half-plane (if not -> not leaving phi extent)
                                                   >> 1803 
                                                   >> 1804       if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  sphi = kInfinity ;
                                                   >> 1805       else     // Leaving via Starting phi
                                                   >> 1806       {
                                                   >> 1807               sidephi = kSPhi  ; 
                                                   >> 1808   
                                                   >> 1809         if ( pDistS > -kCarTolerance*0.5 )  sphi = 0.0 ;
                                                   >> 1810       }
                                                   >> 1811     }
                                                   >> 1812     else sphi = kInfinity ;
                                                   >> 1813   }
                                                   >> 1814   else  // fDPhi <= pi
                                                   >> 1815   {
                                                   >> 1816     if ( compE >= 0.0 )
                                                   >> 1817     {
                                                   >> 1818       if ( compS < 0.0 )
                                                   >> 1819       {
                                                   >> 1820                 
                                                   >> 1821         sphi = pDistS/compS ;
                                                   >> 1822         xi   = p.x() + sphi*v.x() ;
                                                   >> 1823         yi   = p.y() + sphi*v.y() ;
                                                   >> 1824 
                                                   >> 1825 // Check intersection in correct half-plane (if not -> remain in extent)
                                                   >> 1826 
                                                   >> 1827         if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  sphi = kInfinity ;
                                                   >> 1828         else // otherwise leaving via Starting phi
                                                   >> 1829         {
                                                   >> 1830           sidephi = kSPhi ;
                                                   >> 1831         }
                                                   >> 1832       }
                                                   >> 1833       else sphi = kInfinity ;
                                                   >> 1834     }
                                                   >> 1835     else // CompE < 0, leaving immediately by ending
                                                   >> 1836     {
                                                   >> 1837       sidephi = kEPhi ;
                                                   >> 1838       sphi    = 0.0 ;
                                                   >> 1839     }
                                                   >> 1840   }
                                                   >> 1841       }    
1891     }                                            1842     }
1892     else                                         1843     else
1893     {                                            1844     {
1894       // On z axis + travel not || to z axis  << 1845 // On z axis + travel not || to z axis -> if phi of vector direction
1895       // within phi of shape, Step limited by << 1846 // within phi of shape, Step limited by rmax, else Step =0
1896                                                  1847 
1897       if ( (fSPhi-halfAngTolerance <= vphi)   << 1848       vphi = atan2(v.y(),v.x()) ;
1898         && (vphi <= fSPhi+fDPhi+halfAngTolera << 1849 
1899       {                                       << 1850       if ( fSPhi < vphi && vphi < fSPhi + fDPhi ) sphi = kInfinity ;
1900         sphi = kInfinity ;                    << 
1901       }                                       << 
1902       else                                       1851       else
1903       {                                          1852       {
1904         sidephi = kSPhi  ;   // arbitrary        1853         sidephi = kSPhi  ;   // arbitrary 
1905         sphi    = 0.0 ;                       << 1854   sphi    = 0.0 ;
1906       }                                          1855       }
1907     }                                         << 1856     }     
1908     if ( sphi < snxt )  // Order intersecttio    1857     if ( sphi < snxt )  // Order intersecttions
1909     {                                            1858     {
1910       snxt = sphi ;                           << 1859         snxt=sphi ;
1911       side = sidephi ;                        << 1860         side=sidephi ;
1912     }                                            1861     }
1913   }                                              1862   }
1914   if ( srd < snxt )  // Order intersections   << 1863   if ( sr < snxt )  // Order intersections
1915   {                                              1864   {
1916     snxt = srd   ;                            << 1865     snxt = sr    ;
1917     side = sider ;                               1866     side = sider ;
1918   }                                              1867   }
1919   if (calcNorm)                                  1868   if (calcNorm)
1920   {                                              1869   {
1921     switch(side)                                 1870     switch(side)
1922     {                     // Note: returned v << 1871     {
1923       case kRMax:         // (divide by frmax << 1872   case kRMax:
1924         xi         = p.x() + snxt*v.x() ;     << 1873           // Note: returned vector not normalised
1925         yi         = p.y() + snxt*v.y() ;     << 1874           // (divide by frmax for unit vector)
1926         risec      = std::sqrt(xi*xi + yi*yi) << 1875     xi         = p.x() + snxt*v.x() ;
1927         *n         = G4ThreeVector(xi/risec,y << 1876     yi         = p.y() + snxt*v.y() ;
1928         *validNorm = true ;                   << 1877     risec      = sqrt(xi*xi + yi*yi)*secRMax ;
1929         break ;                               << 1878     *n         = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1930       case kRMin:                             << 1879     *validNorm = true ;
1931         *validNorm = false ;  // Rmin is inco << 1880     break ;
1932         break ;                               << 1881   case kRMin:
1933       case kSPhi:                             << 1882     *validNorm=false ;  // Rmin is inconvex
1934         if ( fDPhi <= pi )                    << 1883     break ;
1935         {                                     << 1884   case kSPhi:
1936           *n         = G4ThreeVector(sinSPhi, << 1885     if ( fDPhi <= M_PI )
1937           *validNorm = true ;                 << 1886     {
1938         }                                     << 1887       *n         = G4ThreeVector(sin(fSPhi),-cos(fSPhi),0) ;
1939         else                                  << 1888       *validNorm = true ;
1940         {                                     << 1889     }
1941           *validNorm = false ;                << 1890     else   *validNorm = false ;
1942         }                                     << 1891     break ;
1943         break ;                               << 1892   case kEPhi:
1944       case kEPhi:                             << 1893     if ( fDPhi <= M_PI )
1945         if ( fDPhi <= pi )                    << 1894     {
1946         {                                     << 1895       *n         = G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0) ;
1947           *n = G4ThreeVector(-sinEPhi, cosEPh << 1896       *validNorm = true ;
1948           *validNorm = true ;                 << 1897     }
1949         }                                     << 1898     else  *validNorm = false ;
1950         else                                  << 1899     break ;
1951         {                                     << 1900   case kPZ:
1952           *validNorm = false ;                << 1901     *n         = G4ThreeVector(0,0,1) ;
1953         }                                     << 1902     *validNorm = true ;
1954         break ;                               << 1903     break ;
1955       case kPZ:                               << 1904         case kMZ:
1956         *n         = G4ThreeVector(0,0,1) ;   << 1905     *n         = G4ThreeVector(0,0,-1) ;
1957         *validNorm = true ;                   << 1906     *validNorm = true ;
1958         break ;                               << 1907     break ;
1959       case kMZ:                               << 1908   default:
1960         *n         = G4ThreeVector(0,0,-1) ;  << 1909           G4cout.precision(16) ;
1961         *validNorm = true ;                   << 1910           G4cout << G4endl ;
1962         break ;                               << 1911           G4cout << "Cons parameters:" << G4endl << G4endl ;
1963       default:                                << 1912           G4cout << "fRmin1 = "  << fRmin1/mm << " mm" << G4endl ;
1964         G4cout << G4endl ;                    << 1913           G4cout << "fRmax1 = "  << fRmax1/mm << " mm" << G4endl ;
1965         DumpInfo();                           << 1914           G4cout << "fRmin2 = "  << fRmin2/mm << " mm" << G4endl ;
1966         std::ostringstream message;           << 1915           G4cout << "fRmax2 = "  << fRmax2/mm << " mm" << G4endl ;
1967         G4long oldprc = message.precision(16) << 1916           G4cout << "fDz = "     << fDz/mm << " mm" << G4endl ;
1968         message << "Undefined side for valid  << 1917           G4cout << "fSPhi = "   << fSPhi/degree << " degree" << G4endl ;
1969                 << G4endl                     << 1918           G4cout << "fDPhi = "   << fDPhi/degree << " degree" << G4endl ;
1970                 << "Position:"  << G4endl <<  << 1919           G4cout << "Position:"  << G4endl << G4endl ;
1971                 << "p.x() = "   << p.x()/mm < << 1920           G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1972                 << "p.y() = "   << p.y()/mm < << 1921           G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1973                 << "p.z() = "   << p.z()/mm < << 1922           G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1974                 << "pho at z = "   << std::sq << 1923           G4cout << "Direction:" << G4endl << G4endl ;
1975                 << " mm" << G4endl << G4endl  << 1924           G4cout << "v.x() = "   << v.x() << G4endl ;
1976         if( p.x() != 0. || p.y() != 0.)       << 1925           G4cout << "v.y() = "   << v.y() << G4endl ;
1977         {                                     << 1926           G4cout << "v.z() = "   << v.z() << G4endl<< G4endl ;
1978            message << "point phi = "   << std << 1927           G4cout << "Proposed distance :" << G4endl<< G4endl ;
1979                    << " degrees" << G4endl << << 1928           G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl ;
1980         }                                     << 1929           G4Exception("Invalid enum in G4Cons::DistanceToOut") ;
1981         message << "Direction:" << G4endl <<  << 1930     break ;
1982                 << "v.x() = "   << v.x() << G << 
1983                 << "v.y() = "   << v.y() << G << 
1984                 << "v.z() = "   << v.z() << G << 
1985                 << "Proposed distance :" << G << 
1986                 << "snxt = "    << snxt/mm << << 
1987         message.precision(oldprc) ;           << 
1988         G4Exception("G4Cons::DistanceToOut(p, << 
1989                     JustWarning, message) ;   << 
1990         break ;                               << 
1991     }                                            1931     }
1992   }                                              1932   }
1993   if (snxt < halfCarTolerance)  { snxt = 0.;  << 1933   if (snxt < kCarTolerance*0.5) snxt = 0.;
1994                                               << 
1995   return snxt ;                                  1934   return snxt ;
1996 }                                                1935 }
1997                                                  1936 
1998 /////////////////////////////////////////////    1937 //////////////////////////////////////////////////////////////////
1999 //                                               1938 //
2000 // Calculate distance (<=actual) to closest s    1939 // Calculate distance (<=actual) to closest surface of shape from inside
2001                                                  1940 
2002 G4double G4Cons::DistanceToOut(const G4ThreeV    1941 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
2003 {                                                1942 {
2004   G4double safe=0.0, rho, safeR1, safeR2, saf << 1943   G4double safe,rho,safeR1,safeR2,safeZ ;
2005   G4double tanRMin, secRMin, pRMin;           << 1944   G4double tanRMin,secRMin,pRMin ;
2006   G4double tanRMax, secRMax, pRMax;           << 1945   G4double tanRMax,secRMax,pRMax ;
                                                   >> 1946   G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi ;
2007                                                  1947 
2008 #ifdef G4CSGDEBUG                                1948 #ifdef G4CSGDEBUG
2009   if( Inside(p) == kOutside )                    1949   if( Inside(p) == kOutside )
2010   {                                              1950   {
2011     G4int oldprc=G4cout.precision(16) ;       << 1951     G4cout.precision(16) ;
2012     G4cout << G4endl ;                           1952     G4cout << G4endl ;
2013     DumpInfo();                               << 1953     G4cout << "Cons parameters:" << G4endl << G4endl ;
                                                   >> 1954     G4cout << "fRmin1 = "  << fRmin1/mm << " mm" << G4endl ;
                                                   >> 1955     G4cout << "fRmax1 = "  << fRmax1/mm << " mm" << G4endl ;
                                                   >> 1956     G4cout << "fRmin2 = "  << fRmin2/mm << " mm" << G4endl ;
                                                   >> 1957     G4cout << "fRmax2 = "  << fRmax2/mm << " mm" << G4endl ;
                                                   >> 1958     G4cout << "fDz = "     << fDz/mm << " mm" << G4endl ;
                                                   >> 1959     G4cout << "fSPhi = "   << fSPhi/degree << " degree" << G4endl ;
                                                   >> 1960     G4cout << "fDPhi = "   << fDPhi/degree << " degree" << G4endl << G4endl ;
2014     G4cout << "Position:"  << G4endl << G4end    1961     G4cout << "Position:"  << G4endl << G4endl ;
2015     G4cout << "p.x() = "   << p.x()/mm << " m    1962     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
2016     G4cout << "p.y() = "   << p.y()/mm << " m    1963     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
2017     G4cout << "p.z() = "   << p.z()/mm << " m    1964     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
2018     G4cout << "pho at z = "   << std::sqrt( p << 1965     G4cout << "G4Cons::DistanceToOut(p) - point p is outside ?!" << G4endl ;
2019            << " mm" << G4endl << G4endl ;     << 1966     // G4Exception("Invalid call in G4Cons::DistanceToOut(p), p is outside") ;
2020     if( (p.x() != 0.) || (p.x() != 0.) )      << 
2021     {                                         << 
2022       G4cout << "point phi = "   << std::atan << 
2023              << " degrees" << G4endl << G4end << 
2024     }                                         << 
2025     G4cout.precision(oldprc) ;                << 
2026     G4Exception("G4Cons::DistanceToOut(p)", " << 
2027                 JustWarning, "Point p is outs << 
2028   }                                              1967   }
2029 #endif                                           1968 #endif
2030                                                  1969 
2031   rho = std::sqrt(p.x()*p.x() + p.y()*p.y())  << 1970   rho = sqrt(p.x()*p.x() + p.y()*p.y()) ;
2032   safeZ = fDz - std::fabs(p.z()) ;            << 1971   safeZ = fDz - fabs(p.z()) ;
2033                                                  1972 
2034   if ((fRmin1 != 0.0) || (fRmin2 != 0.0))     << 1973   if (fRmin1 || fRmin2)
2035   {                                              1974   {
2036     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;        1975     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2037     secRMin = std::sqrt(1.0 + tanRMin*tanRMin << 1976     secRMin = sqrt(1.0 + tanRMin*tanRMin) ;
2038     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin    1977     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2039     safeR1  = (rho - pRMin)/secRMin ;            1978     safeR1  = (rho - pRMin)/secRMin ;
2040   }                                              1979   }
2041   else                                        << 1980   else safeR1 = kInfinity ;
2042   {                                           << 
2043     safeR1 = kInfinity ;                      << 
2044   }                                           << 
2045                                                  1981 
2046   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;          1982   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2047   secRMax = std::sqrt(1.0 + tanRMax*tanRMax)  << 1983   secRMax = sqrt(1.0 + tanRMax*tanRMax) ;
2048   pRMax   = tanRMax*p.z() + (fRmax1+fRmax2)*0    1984   pRMax   = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2049   safeR2  = (pRMax - rho)/secRMax ;              1985   safeR2  = (pRMax - rho)/secRMax ;
2050                                                  1986 
2051   if (safeR1 < safeR2)  { safe = safeR1; }    << 1987   if (safeR1 < safeR2) safe = safeR1 ;
2052   else                  { safe = safeR2; }    << 1988   else                 safe = safeR2 ;
2053   if (safeZ < safe)     { safe = safeZ ; }    << 1989   if (safeZ < safe)    safe = safeZ  ;
2054                                                  1990 
2055   // Check if phi divided, Calc distances clo << 1991 // Check if phi divided, Calc distances closest phi plane
2056                                                  1992 
2057   if (!fPhiFullCone)                          << 1993   if (fDPhi < 2.0*M_PI)
2058   {                                              1994   {
2059     // Above/below central phi of G4Cons?     << 1995 // Above/below central phi of G4Cons?
                                                   >> 1996 
                                                   >> 1997     phiC    = fSPhi + fDPhi*0.5 ;
                                                   >> 1998     cosPhiC = cos(phiC) ;
                                                   >> 1999     sinPhiC = sin(phiC) ;
2060                                                  2000 
2061     if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 << 2001     if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
2062     {                                            2002     {
2063       safePhi = -(p.x()*sinSPhi - p.y()*cosSP << 2003       safePhi = -(p.x()*sin(fSPhi) - p.y()*cos(fSPhi)) ;
2064     }                                            2004     }
2065     else                                         2005     else
2066     {                                            2006     {
2067       safePhi = (p.x()*sinEPhi - p.y()*cosEPh << 2007       ePhi    = fSPhi + fDPhi ;
                                                   >> 2008       safePhi = (p.x()*sin(ePhi) - p.y()*cos(ePhi)) ;
2068     }                                            2009     }
2069     if (safePhi < safe)  { safe = safePhi; }  << 2010     if (safePhi < safe) safe = safePhi ;
2070   }                                              2011   }
2071   if ( safe < 0 )  { safe = 0; }              << 2012   if ( safe < 0 ) safe = 0 ;
2072                                               << 2013   return safe ; 
2073   return safe ;                               << 
2074 }                                                2014 }
2075                                                  2015 
2076 ///////////////////////////////////////////// << 2016 ////////////////////////////////////////////////////////////////////////////
2077 //                                            << 
2078 // GetEntityType                              << 
2079                                               << 
2080 G4GeometryType G4Cons::GetEntityType() const  << 
2081 {                                             << 
2082   return {"G4Cons"};                          << 
2083 }                                             << 
2084                                               << 
2085 ///////////////////////////////////////////// << 
2086 //                                            << 
2087 // Make a clone of the object                 << 
2088 //                                            << 
2089 G4VSolid* G4Cons::Clone() const               << 
2090 {                                             << 
2091   return new G4Cons(*this);                   << 
2092 }                                             << 
2093                                               << 
2094 ///////////////////////////////////////////// << 
2095 //                                               2017 //
2096 // Stream object contents to an output stream << 2018 // Create a List containing the transformed vertices
                                                   >> 2019 // Ordering [0-3] -fDz cross section
                                                   >> 2020 //          [4-7] +fDz cross section such that [0] is below [4],
                                                   >> 2021 //                                             [1] below [5] etc.
                                                   >> 2022 // Note:
                                                   >> 2023 //  Caller has deletion resposibility
                                                   >> 2024 //  Potential improvement: For last slice, use actual ending angle
                                                   >> 2025 //                         to avoid rounding error problems.
2097                                                  2026 
2098 std::ostream& G4Cons::StreamInfo(std::ostream << 2027 G4ThreeVectorList*
                                                   >> 2028 G4Cons::CreateRotatedVertices(const G4AffineTransform& pTransform) const
2099 {                                                2029 {
2100   G4long oldprc = os.precision(16);           << 2030   G4ThreeVectorList* vertices ;
2101   os << "------------------------------------ << 2031   G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
2102      << "    *** Dump for solid - " << GetNam << 2032   G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
2103      << "    ================================ << 2033   G4double cosCrossAngle, sinCrossAngle, sAngle ;
2104      << " Solid type: G4Cons\n"               << 2034   G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
2105      << " Parameters: \n"                     << 2035   G4int crossSection, noCrossSections ;
2106      << "   inside  -fDz radius: "  << fRmin1 << 
2107      << "   outside -fDz radius: "  << fRmax1 << 
2108      << "   inside  +fDz radius: "  << fRmin2 << 
2109      << "   outside +fDz radius: "  << fRmax2 << 
2110      << "   half length in Z   : "  << fDz/mm << 
2111      << "   starting angle of segment: " << f << 
2112      << "   delta angle of segment   : " << f << 
2113      << "------------------------------------ << 
2114   os.precision(oldprc);                       << 
2115                                                  2036 
2116   return os;                                  << 2037 // Compute no of cross-sections necessary to mesh cone
2117 }                                             << 2038     
                                                   >> 2039   noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
2118                                                  2040 
                                                   >> 2041   if (noCrossSections < kMinMeshSections)
                                                   >> 2042   {
                                                   >> 2043     noCrossSections = kMinMeshSections ;
                                                   >> 2044   }
                                                   >> 2045   else if (noCrossSections > kMaxMeshSections)
                                                   >> 2046   {
                                                   >> 2047     noCrossSections = kMaxMeshSections ;
                                                   >> 2048   }
                                                   >> 2049   meshAngle = fDPhi/(noCrossSections - 1) ;
2119                                                  2050 
                                                   >> 2051 // G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1  ;
2120                                                  2052 
2121 ///////////////////////////////////////////// << 2053   meshRMax1 = fRmax1/cos(meshAngle*0.5) ;
2122 //                                            << 2054   meshRMax2 = fRmax2/cos(meshAngle*0.5) ;
2123 // GetPointOnSurface                          << 
2124                                                  2055 
2125 G4ThreeVector G4Cons::GetPointOnSurface() con << 2056 // If complete in phi, set start angle such that mesh will be at RMax
2126 {                                             << 2057 // on the x axis. Will give better extent calculations when not rotated.
2127   // declare working variables                << 
2128   //                                          << 
2129   G4double rone = (fRmax1-fRmax2)/(2.*fDz);   << 
2130   G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);   << 
2131   G4double qone = (fRmax1 == fRmax2) ? 0. : f << 
2132   G4double qtwo = (fRmin1 == fRmin2) ? 0. : f << 
2133                                               << 
2134   G4double slin   = std::hypot(fRmin1-fRmin2, << 
2135   G4double slout  = std::hypot(fRmax1-fRmax2, << 
2136   G4double Aone   = 0.5*fDPhi*(fRmax2 + fRmax << 
2137   G4double Atwo   = 0.5*fDPhi*(fRmin2 + fRmin << 
2138   G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1- << 
2139   G4double Afour  = 0.5*fDPhi*(fRmax2*fRmax2- << 
2140   G4double Afive  = fDz*(fRmax1-fRmin1+fRmax2 << 
2141                                               << 
2142   G4double phi    = G4RandFlat::shoot(fSPhi,f << 
2143   G4double cosu   = std::cos(phi);            << 
2144   G4double sinu   = std::sin(phi);            << 
2145   G4double rRand1 = GetRadiusInRing(fRmin1, f << 
2146   G4double rRand2 = GetRadiusInRing(fRmin2, f << 
2147                                               << 
2148   if ( (fSPhi == 0.) && fPhiFullCone )  { Afi << 
2149   G4double chose  = G4RandFlat::shoot(0.,Aone << 
2150                                                  2058 
2151   if( (chose >= 0.) && (chose < Aone) ) // ou << 2059   if (fDPhi == M_PI*2.0 && fSPhi == 0.0 )
2152   {                                              2060   {
2153     if(fRmax1 != fRmax2)                      << 2061     sAngle = -meshAngle*0.5 ;
2154     {                                         << 
2155       G4double zRand = G4RandFlat::shoot(-1.* << 
2156       return { rone*cosu*(qone-zRand), rone*s << 
2157     }                                         << 
2158     else                                      << 
2159     {                                         << 
2160       return { fRmax1*cosu, fRmax2*sinu, G4Ra << 
2161     }                                         << 
2162   }                                              2062   }
2163   else if( (chose >= Aone) && (chose < Aone + << 2063   else
2164   {                                              2064   {
2165     if(fRmin1 != fRmin2)                      << 2065     sAngle = fSPhi ;
2166     {                                         << 2066   } 
2167       G4double zRand = G4RandFlat::shoot(-1.* << 2067   vertices = new G4ThreeVectorList();
2168       return { rtwo*cosu*(qtwo-zRand), rtwo*s << 2068   vertices->reserve(noCrossSections*4) ;
2169     }                                         << 2069 
2170     else                                      << 2070   if (vertices)
2171     {                                         << 2071   {
2172       return { fRmin1*cosu, fRmin2*sinu, G4Ra << 2072     for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
                                                   >> 2073     {
                                                   >> 2074 // Compute coordinates of cross section at section crossSection
                                                   >> 2075 
                                                   >> 2076       crossAngle    = sAngle + crossSection*meshAngle ;
                                                   >> 2077       cosCrossAngle = cos(crossAngle) ;
                                                   >> 2078       sinCrossAngle = sin(crossAngle) ;
                                                   >> 2079 
                                                   >> 2080       rMaxX1 = meshRMax1*cosCrossAngle ;
                                                   >> 2081       rMaxY1 = meshRMax1*sinCrossAngle ;
                                                   >> 2082       rMaxX2 = meshRMax2*cosCrossAngle ;
                                                   >> 2083       rMaxY2 = meshRMax2*sinCrossAngle ;
                                                   >> 2084         
                                                   >> 2085 // G4double RMin = (fRmin2 <= fRmin1) ? fRmin2 : fRmin1  ;
                                                   >> 2086 
                                                   >> 2087       rMinX1 = fRmin1*cosCrossAngle ;
                                                   >> 2088       rMinY1 = fRmin1*sinCrossAngle ;
                                                   >> 2089       rMinX2 = fRmin2*cosCrossAngle ;
                                                   >> 2090       rMinY2 = fRmin2*sinCrossAngle ;
                                                   >> 2091         
                                                   >> 2092       vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
                                                   >> 2093       vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
                                                   >> 2094       vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
                                                   >> 2095       vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;
                                                   >> 2096 
                                                   >> 2097       vertices->push_back(pTransform.TransformPoint(vertex0)) ;
                                                   >> 2098       vertices->push_back(pTransform.TransformPoint(vertex1)) ;
                                                   >> 2099       vertices->push_back(pTransform.TransformPoint(vertex2)) ;
                                                   >> 2100       vertices->push_back(pTransform.TransformPoint(vertex3)) ;
2173     }                                            2101     }
2174   }                                              2102   }
2175   else if( (chose >= Aone + Atwo) && (chose < << 2103   else
2176   {                                           << 
2177     return {rRand1*cosu, rRand1*sinu, -1*fDz} << 
2178   }                                           << 
2179   else if( (chose >= Aone + Atwo + Athree)    << 
2180         && (chose < Aone + Atwo + Athree + Af << 
2181   {                                           << 
2182     return { rRand2*cosu, rRand2*sinu, fDz }; << 
2183   }                                           << 
2184   else if( (chose >= Aone + Atwo + Athree + A << 
2185         && (chose < Aone + Atwo + Athree + Af << 
2186   {                                           << 
2187     G4double zRand  = G4RandFlat::shoot(-1.*f << 
2188     rRand1 = G4RandFlat::shoot(fRmin2-((zRand << 
2189                                fRmax2-((zRand << 
2190     return { rRand1*cosSPhi, rRand1*sinSPhi,  << 
2191   }                                           << 
2192   else // SPhi+DPhi section                   << 
2193   {                                              2104   {
2194     G4double zRand  = G4RandFlat::shoot(-1.*f << 2105     G4Exception("G4Cons::CreateRotatedVertices - Out of memory !");
2195     rRand1 = G4RandFlat::shoot(fRmin2-((zRand << 
2196                                fRmax2-((zRand << 
2197     return { rRand1*cosEPhi, rRand1*sinEPhi,  << 
2198   }                                              2106   }
                                                   >> 2107   return vertices ;
2199 }                                                2108 }
2200                                                  2109 
2201 /////////////////////////////////////////////    2110 //////////////////////////////////////////////////////////////////////////
2202 //                                               2111 //
2203 // Methods for visualisation                     2112 // Methods for visualisation
2204                                                  2113 
2205 void G4Cons::DescribeYourselfTo (G4VGraphicsS    2114 void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const
2206 {                                                2115 {
2207   scene.AddSolid (*this);                     << 2116   scene.AddThis (*this);
2208 }                                                2117 }
2209                                                  2118 
2210 G4Polyhedron* G4Cons::CreatePolyhedron () con    2119 G4Polyhedron* G4Cons::CreatePolyhedron () const
2211 {                                                2120 {
2212   return new G4PolyhedronCons(fRmin1,fRmax1,f << 2121     return new G4PolyhedronCons(fRmin1, fRmax1, fRmin2, fRmax2, fDz, fSPhi,
                                                   >> 2122                                 fDPhi);
2213 }                                                2123 }
2214                                                  2124 
2215 #endif                                        << 2125 G4NURBS* G4Cons::CreateNURBS () const
                                                   >> 2126 {
                                                   >> 2127     G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
                                                   >> 2128     
                                                   >> 2129     return new G4NURBSbox (RMax, RMax, fDz);       // Box for now!!!
                                                   >> 2130 }
                                                   >> 2131 
                                                   >> 2132 //
                                                   >> 2133 //
                                                   >> 2134 /////////////////////////////// End of G4Cons.cc file //////////////////////
2216                                                  2135