source: trunk/source/geometry/solids/CSG/src/G4Tubs.cc @ 1358

Last change on this file since 1358 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 54.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4Tubs.cc,v 1.79 2009/06/30 10:10:11 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30//
31// class G4Tubs
32//
33// History:
34//
35// 02.08.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for negative value under sqrt
36//                      for the case: p on the surface and v is tangent to the surface
37// 11.05.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for phi < 2pi
38// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
39// 16.03.05 V.Grichine: SurfaceNormal(p) with edges/corners for boolean
40// 20.07.01 V.Grichine: bug fixed in Inside(p)
41// 20.02.01 V.Grichine: bug fixed in Inside(p) and CalculateExtent was
42//                      simplified base on G4Box::CalculateExtent
43// 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
44// 28.11.00 V.Grichine: bug fixed in Inside(p)
45// 31.10.00 V.Grichine: assign sr, sphi in Distance ToOut(p,v,...)
46// 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
47// 02.08.00 V.Grichine: point is outside check in Distance ToOut(p)
48// 17.05.00 V.Grichine: bugs (#76,#91) fixed in Distance ToOut(p,v,...)
49// 31.03.00 V.Grichine: bug fixed in Inside(p)
50// 19.11.99 V.Grichine: side = kNull in DistanceToOut(p,v,...)
51// 13.10.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
52// 28.05.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
53// 25.05.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
54// 23.03.99 V.Grichine: bug fixed in DistanceToIn(p,v)
55// 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
56// 18.06.98 V.Grichine: n-normalisation in DistanceToOut(p,v)
57//
58// 1994-95  P.Kent:     implementation
59//
60/////////////////////////////////////////////////////////////////////////
61
62#include "G4Tubs.hh"
63
64#include "G4VoxelLimits.hh"
65#include "G4AffineTransform.hh"
66#include "G4GeometryTolerance.hh"
67
68#include "G4VPVParameterisation.hh"
69
70#include "Randomize.hh"
71
72#include "meshdefs.hh"
73
74#include "G4VGraphicsScene.hh"
75#include "G4Polyhedron.hh"
76#include "G4NURBS.hh"
77#include "G4NURBStube.hh"
78#include "G4NURBScylinder.hh"
79#include "G4NURBStubesector.hh"
80
81using namespace CLHEP;
82
83/////////////////////////////////////////////////////////////////////////
84//
85// Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
86//             - note if pdphi>2PI then reset to 2PI
87
88G4Tubs::G4Tubs( const G4String &pName,
89                      G4double pRMin, G4double pRMax,
90                      G4double pDz,
91                      G4double pSPhi, G4double pDPhi )
92  : G4CSGSolid(pName), fSPhi(0), fDPhi(0)
93{
94
95  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
96  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
97
98  if (pDz>0) // Check z-len
99  {
100    fDz = pDz ;
101  }
102  else
103  {
104    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
105           << "        Negative Z half-length (" << pDz << ") in solid: "
106           << GetName() << G4endl;
107    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
108                "Invalid Z half-length");
109  }
110  if ( (pRMin < pRMax) && (pRMin >= 0) ) // Check radii
111  {
112    fRMin = pRMin ; 
113    fRMax = pRMax ;
114  }
115  else
116  {
117    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
118           << "        Invalid values for radii in solid " << GetName()
119           << G4endl
120           << "        pRMin = " << pRMin << ", pRMax = " << pRMax << G4endl;
121    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
122                "Invalid radii.");
123  }
124
125  // Check angles
126
127  CheckPhiAngles(pSPhi, pDPhi);
128}
129
130///////////////////////////////////////////////////////////////////////
131//
132// Fake default constructor - sets only member data and allocates memory
133//                            for usage restricted to object persistency.
134//
135G4Tubs::G4Tubs( __void__& a )
136  : G4CSGSolid(a)
137{
138}
139
140//////////////////////////////////////////////////////////////////////////
141//
142// Destructor
143
144G4Tubs::~G4Tubs()
145{
146}
147
148/////////////////////////////////////////////////////////////////////////
149//
150// Dispatch to parameterisation for replication mechanism dimension
151// computation & modification.
152
153void G4Tubs::ComputeDimensions(       G4VPVParameterisation* p,
154                                const G4int n,
155                                const G4VPhysicalVolume* pRep )
156{
157  p->ComputeDimensions(*this,n,pRep) ;
158}
159
160////////////////////////////////////////////////////////////////////////
161//
162// Calculate extent under transform and specified limit
163
164G4bool G4Tubs::CalculateExtent( const EAxis              pAxis,
165                                const G4VoxelLimits&     pVoxelLimit,
166                                const G4AffineTransform& pTransform,
167                                      G4double&          pMin, 
168                                      G4double&          pMax    ) const
169{
170
171  if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
172  {
173    // Special case handling for unrotated solid tubes
174    // Compute x/y/z mins and maxs fro bounding box respecting limits,
175    // with early returns if outside limits. Then switch() on pAxis,
176    // and compute exact x and y limit for x/y case
177     
178    G4double xoffset, xMin, xMax;
179    G4double yoffset, yMin, yMax;
180    G4double zoffset, zMin, zMax;
181
182    G4double diff1, diff2, maxDiff, newMin, newMax;
183    G4double xoff1, xoff2, yoff1, yoff2, delta;
184
185    xoffset = pTransform.NetTranslation().x();
186    xMin = xoffset - fRMax;
187    xMax = xoffset + fRMax;
188
189    if (pVoxelLimit.IsXLimited())
190    {
191      if ( (xMin > pVoxelLimit.GetMaxXExtent())
192        || (xMax < pVoxelLimit.GetMinXExtent()) )
193      {
194        return false;
195      }
196      else
197      {
198        if (xMin < pVoxelLimit.GetMinXExtent())
199        {
200          xMin = pVoxelLimit.GetMinXExtent();
201        }
202        if (xMax > pVoxelLimit.GetMaxXExtent())
203        {
204          xMax = pVoxelLimit.GetMaxXExtent();
205        }
206      }
207    }
208    yoffset = pTransform.NetTranslation().y();
209    yMin    = yoffset - fRMax;
210    yMax    = yoffset + fRMax;
211
212    if ( pVoxelLimit.IsYLimited() )
213    {
214      if ( (yMin > pVoxelLimit.GetMaxYExtent())
215        || (yMax < pVoxelLimit.GetMinYExtent()) )
216      {
217        return false;
218      }
219      else
220      {
221        if (yMin < pVoxelLimit.GetMinYExtent())
222        {
223          yMin = pVoxelLimit.GetMinYExtent();
224        }
225        if (yMax > pVoxelLimit.GetMaxYExtent())
226        {
227          yMax=pVoxelLimit.GetMaxYExtent();
228        }
229      }
230    }
231    zoffset = pTransform.NetTranslation().z();
232    zMin    = zoffset - fDz;
233    zMax    = zoffset + fDz;
234
235    if ( pVoxelLimit.IsZLimited() )
236    {
237      if ( (zMin > pVoxelLimit.GetMaxZExtent())
238        || (zMax < pVoxelLimit.GetMinZExtent()) )
239      {
240        return false;
241      }
242      else
243      {
244        if (zMin < pVoxelLimit.GetMinZExtent())
245        {
246          zMin = pVoxelLimit.GetMinZExtent();
247        }
248        if (zMax > pVoxelLimit.GetMaxZExtent())
249        {
250          zMax = pVoxelLimit.GetMaxZExtent();
251        }
252      }
253    }
254    switch ( pAxis )  // Known to cut cylinder
255    {
256      case kXAxis :
257      {
258        yoff1 = yoffset - yMin;
259        yoff2 = yMax    - yoffset;
260
261        if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
262        {                                   // => no change
263          pMin = xMin;
264          pMax = xMax;
265        }
266        else
267        {
268          // Y limits don't cross max/min x => compute max delta x,
269          // hence new mins/maxs
270
271          delta   = fRMax*fRMax - yoff1*yoff1;
272          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
273          delta   = fRMax*fRMax - yoff2*yoff2;
274          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
275          maxDiff = (diff1 > diff2) ? diff1:diff2;
276          newMin  = xoffset - maxDiff;
277          newMax  = xoffset + maxDiff;
278          pMin    = (newMin < xMin) ? xMin : newMin;
279          pMax    = (newMax > xMax) ? xMax : newMax;
280        }   
281        break;
282      }
283      case kYAxis :
284      {
285        xoff1 = xoffset - xMin;
286        xoff2 = xMax - xoffset;
287
288        if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
289        {                                   // => no change
290          pMin = yMin;
291          pMax = yMax;
292        }
293        else
294        {
295          // X limits don't cross max/min y => compute max delta y,
296          // hence new mins/maxs
297
298          delta   = fRMax*fRMax - xoff1*xoff1;
299          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
300          delta   = fRMax*fRMax - xoff2*xoff2;
301          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
302          maxDiff = (diff1 > diff2) ? diff1 : diff2;
303          newMin  = yoffset - maxDiff;
304          newMax  = yoffset + maxDiff;
305          pMin    = (newMin < yMin) ? yMin : newMin;
306          pMax    = (newMax > yMax) ? yMax : newMax;
307        }
308        break;
309      }
310      case kZAxis:
311      {
312        pMin = zMin;
313        pMax = zMax;
314        break;
315      }
316      default:
317        break;
318    }
319    pMin -= kCarTolerance;
320    pMax += kCarTolerance;
321    return true;
322  }
323  else // Calculate rotated vertex coordinates
324  {
325    G4int i, noEntries, noBetweenSections4;
326    G4bool existsAfterClip = false;
327    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
328
329    pMin =  kInfinity;
330    pMax = -kInfinity;
331
332    noEntries = vertices->size();
333    noBetweenSections4 = noEntries - 4;
334   
335    for ( i = 0 ; i < noEntries ; i += 4 )
336    {
337      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
338    }
339    for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
340    {
341      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
342    }
343    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
344    {
345      existsAfterClip = true;
346      pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
347      pMax += kCarTolerance;
348    }
349    else
350    {
351      // Check for case where completely enveloping clipping volume
352      // If point inside then we are confident that the solid completely
353      // envelopes the clipping volume. Hence set min/max extents according
354      // to clipping volume extents along the specified axis.
355
356      G4ThreeVector clipCentre(
357             (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
358             (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
359             (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
360       
361      if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
362      {
363        existsAfterClip = true;
364        pMin            = pVoxelLimit.GetMinExtent(pAxis);
365        pMax            = pVoxelLimit.GetMaxExtent(pAxis);
366      }
367    }
368    delete vertices;
369    return existsAfterClip;
370  }
371}
372
373
374///////////////////////////////////////////////////////////////////////////
375//
376// Return whether point inside/outside/on surface
377
378EInside G4Tubs::Inside( const G4ThreeVector& p ) const
379{
380  G4double r2,pPhi,tolRMin,tolRMax;
381  EInside in = kOutside ;
382  static const G4double halfCarTolerance=kCarTolerance*0.5;
383  static const G4double halfRadTolerance=kRadTolerance*0.5;
384  static const G4double halfAngTolerance=kAngTolerance*0.5;
385
386  if (std::fabs(p.z()) <= fDz - halfCarTolerance)
387  {
388    r2 = p.x()*p.x() + p.y()*p.y() ;
389
390    if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
391    else       { tolRMin = 0 ; }
392
393    tolRMax = fRMax - halfRadTolerance ;
394     
395    if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
396    {
397      if ( fPhiFullTube )
398      {
399        in = kInside ;
400      }
401      else
402      {
403        // Try inner tolerant phi boundaries (=>inside)
404        // if not inside, try outer tolerant phi boundaries
405
406        if ((tolRMin==0)&&(p.x()<=halfCarTolerance)&&(p.y()<=halfCarTolerance))
407        {
408          in=kSurface;
409        }
410        else
411        {
412          pPhi = std::atan2(p.y(),p.x()) ;
413          if ( pPhi < -halfAngTolerance )  { pPhi += twopi; } // 0<=pPhi<2pi
414
415          if ( fSPhi >= 0 )
416          {
417            if ( (std::abs(pPhi) < halfAngTolerance)
418              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
419            { 
420              pPhi += twopi ; // 0 <= pPhi < 2pi
421            }
422            if ( (pPhi >= fSPhi + halfAngTolerance)
423              && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
424            {
425              in = kInside ;
426            }
427            else if ( (pPhi >= fSPhi - halfAngTolerance)
428                   && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
429            {
430              in = kSurface ;
431            }
432          }
433          else  // fSPhi < 0
434          {
435            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
436              && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;} //kOutside
437            else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
438                   && (pPhi >= fSPhi + fDPhi  - halfAngTolerance) )
439            {
440              in = kSurface ;
441            }
442            else
443            {
444              in = kInside ;
445            }
446          }
447        }                   
448      }
449    }
450    else  // Try generous boundaries
451    {
452      tolRMin = fRMin - halfRadTolerance ;
453      tolRMax = fRMax + halfRadTolerance ;
454
455      if ( tolRMin < 0 )  { tolRMin = 0; }
456
457      if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
458      {
459        if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
460        {                        // Continuous in phi or on z-axis
461          in = kSurface ;
462        }
463        else // Try outer tolerant phi boundaries only
464        {
465          pPhi = std::atan2(p.y(),p.x()) ;
466
467          if ( pPhi < -halfAngTolerance)  { pPhi += twopi; } // 0<=pPhi<2pi
468          if ( fSPhi >= 0 )
469          {
470            if ( (std::abs(pPhi) < halfAngTolerance)
471              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
472            { 
473              pPhi += twopi ; // 0 <= pPhi < 2pi
474            }
475            if ( (pPhi >= fSPhi - halfAngTolerance)
476              && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
477            {
478              in = kSurface ;
479            }
480          }
481          else  // fSPhi < 0
482          {
483            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
484              && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
485            else
486            {
487              in = kSurface ;
488            }
489          }
490        }
491      }
492    }
493  }
494  else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
495  {                                          // Check within tolerant r limits
496    r2      = p.x()*p.x() + p.y()*p.y() ;
497    tolRMin = fRMin - halfRadTolerance ;
498    tolRMax = fRMax + halfRadTolerance ;
499
500    if ( tolRMin < 0 )  { tolRMin = 0; }
501
502    if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
503    {
504      if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
505      {                        // Continuous in phi or on z-axis
506        in = kSurface ;
507      }
508      else // Try outer tolerant phi boundaries
509      {
510        pPhi = std::atan2(p.y(),p.x()) ;
511
512        if ( pPhi < -halfAngTolerance )  { pPhi += twopi; }  // 0<=pPhi<2pi
513        if ( fSPhi >= 0 )
514        {
515          if ( (std::abs(pPhi) < halfAngTolerance)
516            && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
517          { 
518            pPhi += twopi ; // 0 <= pPhi < 2pi
519          }
520          if ( (pPhi >= fSPhi - halfAngTolerance)
521            && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
522          {
523            in = kSurface;
524          }
525        }
526        else  // fSPhi < 0
527        {
528          if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
529            && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;}
530          else
531          {
532            in = kSurface ;
533          }
534        }     
535      }
536    }
537  }
538  return in;
539}
540
541///////////////////////////////////////////////////////////////////////////
542//
543// Return unit normal of surface closest to p
544// - note if point on z axis, ignore phi divided sides
545// - unsafe if point close to z axis a rmin=0 - no explicit checks
546
547G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
548{
549  G4int noSurfaces = 0;
550  G4double rho, pPhi;
551  G4double distZ, distRMin, distRMax;
552  G4double distSPhi = kInfinity, distEPhi = kInfinity;
553
554  static const G4double halfCarTolerance = 0.5*kCarTolerance;
555  static const G4double halfAngTolerance = 0.5*kAngTolerance;
556
557  G4ThreeVector norm, sumnorm(0.,0.,0.);
558  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
559  G4ThreeVector nR, nPs, nPe;
560
561  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
562
563  distRMin = std::fabs(rho - fRMin);
564  distRMax = std::fabs(rho - fRMax);
565  distZ    = std::fabs(std::fabs(p.z()) - fDz);
566
567  if (!fPhiFullTube)    // Protected against (0,0,z)
568  {
569    if ( rho > halfCarTolerance )
570    {
571      pPhi = std::atan2(p.y(),p.x());
572   
573      if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
574      else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
575
576      distSPhi = std::fabs(pPhi - fSPhi);       
577      distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
578    }
579    else if( !fRMin )
580    {
581      distSPhi = 0.; 
582      distEPhi = 0.; 
583    }
584    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
585    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
586  }
587  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
588
589  if( distRMax <= halfCarTolerance )
590  {
591    noSurfaces ++;
592    sumnorm += nR;
593  }
594  if( fRMin && (distRMin <= halfCarTolerance) )
595  {
596    noSurfaces ++;
597    sumnorm -= nR;
598  }
599  if( fDPhi < twopi )   
600  {
601    if (distSPhi <= halfAngTolerance) 
602    {
603      noSurfaces ++;
604      sumnorm += nPs;
605    }
606    if (distEPhi <= halfAngTolerance) 
607    {
608      noSurfaces ++;
609      sumnorm += nPe;
610    }
611  }
612  if (distZ <= halfCarTolerance) 
613  {
614    noSurfaces ++;
615    if ( p.z() >= 0.)  { sumnorm += nZ; }
616    else               { sumnorm -= nZ; }
617  }
618  if ( noSurfaces == 0 )
619  {
620#ifdef G4CSGDEBUG
621    G4Exception("G4Tube::SurfaceNormal(p)", "Notification",
622                JustWarning, "Point p is not on surface !?" );
623    G4cout.precision(20);
624    G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
625          << G4endl << G4endl;
626#endif
627     norm = ApproxSurfaceNormal(p);
628  }
629  else if ( noSurfaces == 1 )  { norm = sumnorm; }
630  else                         { norm = sumnorm.unit(); }
631
632  return norm;
633}
634
635/////////////////////////////////////////////////////////////////////////////
636//
637// Algorithm for SurfaceNormal() following the original specification
638// for points not on the surface
639
640G4ThreeVector G4Tubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
641{ ENorm side ;
642  G4ThreeVector norm ;
643  G4double rho, phi ;
644  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
645
646  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
647
648  distRMin = std::fabs(rho - fRMin) ;
649  distRMax = std::fabs(rho - fRMax) ;
650  distZ    = std::fabs(std::fabs(p.z()) - fDz) ;
651
652  if (distRMin < distRMax) // First minimum
653  {
654    if ( distZ < distRMin )
655    {
656       distMin = distZ ;
657       side    = kNZ ;
658    }
659    else
660    {
661      distMin = distRMin ;
662      side    = kNRMin   ;
663    }
664  }
665  else
666  {
667    if ( distZ < distRMax )
668    {
669      distMin = distZ ;
670      side    = kNZ   ;
671    }
672    else
673    {
674      distMin = distRMax ;
675      side    = kNRMax   ;
676    }
677  }   
678  if (!fPhiFullTube  &&  rho ) // Protected against (0,0,z)
679  {
680    phi = std::atan2(p.y(),p.x()) ;
681
682    if ( phi < 0 )  { phi += twopi; }
683
684    if ( fSPhi < 0 )
685    {
686      distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
687    }
688    else
689    {
690      distSPhi = std::fabs(phi - fSPhi)*rho ;
691    }
692    distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
693                                     
694    if (distSPhi < distEPhi) // Find new minimum
695    {
696      if ( distSPhi < distMin )
697      {
698        side = kNSPhi ;
699      }
700    }
701    else
702    {
703      if ( distEPhi < distMin )
704      {
705        side = kNEPhi ;
706      }
707    }
708  }   
709  switch ( side )
710  {
711    case kNRMin : // Inner radius
712    {                     
713      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
714      break ;
715    }
716    case kNRMax : // Outer radius
717    {                 
718      norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
719      break ;
720    }
721    case kNZ : //    + or - dz
722    {                             
723      if ( p.z() > 0 )  { norm = G4ThreeVector(0,0,1) ; }
724      else              { norm = G4ThreeVector(0,0,-1); }
725      break ;
726    }
727    case kNSPhi:
728    {
729      norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
730      break ;
731    }
732    case kNEPhi:
733    {
734      norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
735      break;
736    }
737    default:
738    {
739      DumpInfo();
740      G4Exception("G4Tubs::ApproxSurfaceNormal()", "Notification", JustWarning,
741                  "Undefined side for valid surface normal to solid.");
742      break ;
743    }   
744  }               
745  return norm;
746}
747
748////////////////////////////////////////////////////////////////////
749//
750//
751// Calculate distance to shape from outside, along normalised vector
752// - return kInfinity if no intersection, or intersection distance <= tolerance
753//
754// - Compute the intersection with the z planes
755//        - if at valid r, phi, return
756//
757// -> If point is outer outer radius, compute intersection with rmax
758//        - if at valid phi,z return
759//
760// -> Compute intersection with inner radius, taking largest +ve root
761//        - if valid (in z,phi), save intersction
762//
763//    -> If phi segmented, compute intersections with phi half planes
764//        - return smallest of valid phi intersections and
765//          inner radius intersection
766//
767// NOTE:
768// - 'if valid' implies tolerant checking of intersection points
769
770G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
771                               const G4ThreeVector& v  ) const
772{
773  G4double snxt = kInfinity ;      // snxt = default return value
774  G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
775  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
776  const G4double dRmax = 100.*fRMax;
777
778  static const G4double halfCarTolerance = 0.5*kCarTolerance;
779  static const G4double halfRadTolerance = 0.5*kRadTolerance;
780
781  // Intersection point variables
782  //
783  G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
784  G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables
785 
786  // Calculate tolerant rmin and rmax
787
788  if (fRMin > kRadTolerance)
789  {
790    tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
791    tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
792  }
793  else
794  {
795    tolORMin2 = 0.0 ;
796    tolIRMin2 = 0.0 ;
797  }
798  tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
799  tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
800
801  // Intersection with Z surfaces
802
803  tolIDz = fDz - halfCarTolerance ;
804  tolODz = fDz + halfCarTolerance ;
805
806  if (std::fabs(p.z()) >= tolIDz)
807  {
808    if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
809    {
810      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;   // Z intersect distance
811
812      if(s < 0.0)  { s = 0.0; }
813
814      xi   = p.x() + s*v.x() ;                // Intersection coords
815      yi   = p.y() + s*v.y() ;
816      rho2 = xi*xi + yi*yi ;
817
818      // Check validity of intersection
819
820      if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
821      {
822        if (!fPhiFullTube && rho2)
823        {
824          // Psi = angle made with central (average) phi of shape
825          //
826          inum   = xi*cosCPhi + yi*sinCPhi ;
827          iden   = std::sqrt(rho2) ;
828          cosPsi = inum/iden ;
829          if (cosPsi >= cosHDPhiIT)  { return s ; }
830        }
831        else
832        {
833          return s ;
834        }
835      }
836    }
837    else
838    {
839      if ( snxt<halfCarTolerance )  { snxt=0; }
840      return snxt ;  // On/outside extent, and heading away
841                     // -> cannot intersect
842    }
843  }
844
845  // -> Can not intersect z surfaces
846  //
847  // Intersection with rmax (possible return) and rmin (must also check phi)
848  //
849  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
850  //
851  // Intersects with x^2+y^2=R^2
852  //
853  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
854  //            t1                t2                t3
855
856  t1 = 1.0 - v.z()*v.z() ;
857  t2 = p.x()*v.x() + p.y()*v.y() ;
858  t3 = p.x()*p.x() + p.y()*p.y() ;
859
860  if ( t1 > 0 )        // Check not || to z axis
861  {
862    b = t2/t1 ;
863    c = t3 - fRMax*fRMax ;
864    if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
865    {
866      // Try outer cylinder intersection
867      //          c=(t3-fRMax*fRMax)/t1;
868
869      c /= t1 ;
870      d = b*b - c ;
871
872      if (d >= 0)  // If real root
873      {
874        s = -b - std::sqrt(d) ;
875        if (s >= 0)  // If 'forwards'
876        {
877          if ( s>dRmax ) // Avoid rounding errors due to precision issues on
878          {              // 64 bits systems. Split long distances and recompute
879            G4double fTerm = s-std::fmod(s,dRmax);
880            s = fTerm + DistanceToIn(p+fTerm*v,v);
881          } 
882          // Check z intersection
883          //
884          zi = p.z() + s*v.z() ;
885          if (std::fabs(zi)<=tolODz)
886          {
887            // Z ok. Check phi intersection if reqd
888            //
889            if (fPhiFullTube)
890            {
891              return s ;
892            }
893            else
894            {
895              xi     = p.x() + s*v.x() ;
896              yi     = p.y() + s*v.y() ;
897              cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
898              if (cosPsi >= cosHDPhiIT)  { return s ; }
899            }
900          }  //  end if std::fabs(zi)
901        }    //  end if (s>=0)
902      }      //  end if (d>=0)
903    }        //  end if (r>=fRMax)
904    else 
905    {
906      // Inside outer radius :
907      // check not inside, and heading through tubs (-> 0 to in)
908
909      if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
910      {
911        // Inside both radii, delta r -ve, inside z extent
912
913        if (!fPhiFullTube)
914        {
915          inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
916          iden   = std::sqrt(t3) ;
917          cosPsi = inum/iden ;
918          if (cosPsi >= cosHDPhiIT)
919          {
920            // In the old version, the small negative tangent for the point
921            // on surface was not taken in account, and returning 0.0 ...
922            // New version: check the tangent for the point on surface and
923            // if no intersection, return kInfinity, if intersection instead
924            // return s.
925            //
926            c = t3-fRMax*fRMax; 
927            if ( c<=0.0 )
928            {
929              return 0.0;
930            }
931            else
932            {
933              c = c/t1 ;
934              d = b*b-c;
935              if ( d>=0.0 )
936              {
937                snxt = c/(-b+std::sqrt(d)); // using safe solution
938                                            // for quadratic equation
939                if ( snxt < halfCarTolerance ) { snxt=0; }
940                return snxt ;
941              }     
942              else
943              {
944                return kInfinity;
945              }
946            }
947          } 
948        }
949        else
950        {   
951          // In the old version, the small negative tangent for the point
952          // on surface was not taken in account, and returning 0.0 ...
953          // New version: check the tangent for the point on surface and
954          // if no intersection, return kInfinity, if intersection instead
955          // return s.
956          //
957          c = t3 - fRMax*fRMax; 
958          if ( c<=0.0 )
959          {
960            return 0.0;
961          }
962          else
963          {
964            c = c/t1 ;
965            d = b*b-c;
966            if ( d>=0.0 )
967            {
968              snxt= c/(-b+std::sqrt(d)); // using safe solution
969                                         // for quadratic equation
970              if ( snxt < halfCarTolerance ) { snxt=0; }
971              return snxt ;
972            }     
973            else
974            {
975              return kInfinity;
976            }
977          }
978        } // end if   (!fPhiFullTube)
979      }   // end if   (t3>tolIRMin2)
980    }     // end if   (Inside Outer Radius)
981    if ( fRMin )    // Try inner cylinder intersection
982    {
983      c = (t3 - fRMin*fRMin)/t1 ;
984      d = b*b - c ;
985      if ( d >= 0.0 )  // If real root
986      {
987        // Always want 2nd root - we are outside and know rmax Hit was bad
988        // - If on surface of rmin also need farthest root
989
990        s = -b + std::sqrt(d) ;
991        if (s >= -halfCarTolerance)  // check forwards
992        {
993          // Check z intersection
994          //
995          if(s < 0.0)  { s = 0.0; }
996          if ( s>dRmax ) // Avoid rounding errors due to precision issues seen
997          {              // 64 bits systems. Split long distances and recompute
998            G4double fTerm = s-std::fmod(s,dRmax);
999            s = fTerm + DistanceToIn(p+fTerm*v,v);
1000          } 
1001          zi = p.z() + s*v.z() ;
1002          if (std::fabs(zi) <= tolODz)
1003          {
1004            // Z ok. Check phi
1005            //
1006            if ( fPhiFullTube )
1007            {
1008              return s ; 
1009            }
1010            else
1011            {
1012              xi     = p.x() + s*v.x() ;
1013              yi     = p.y() + s*v.y() ;
1014              cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1015              if (cosPsi >= cosHDPhiIT)
1016              {
1017                // Good inner radius isect
1018                // - but earlier phi isect still possible
1019
1020                snxt = s ;
1021              }
1022            }
1023          }        //    end if std::fabs(zi)
1024        }          //    end if (s>=0)
1025      }            //    end if (d>=0)
1026    }              //    end if (fRMin)
1027  }
1028
1029  // Phi segment intersection
1030  //
1031  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1032  //
1033  // o NOTE: Large duplication of code between sphi & ephi checks
1034  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1035  //            intersection check <=0 -> >=0
1036  //         -> use some form of loop Construct ?
1037  //
1038  if ( !fPhiFullTube )
1039  {
1040    // First phi surface (Starting phi)
1041    //
1042    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
1043                   
1044    if ( Comp < 0 )  // Component in outwards normal dirn
1045    {
1046      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1047
1048      if ( Dist < halfCarTolerance )
1049      {
1050        s = Dist/Comp ;
1051
1052        if (s < snxt)
1053        {
1054          if ( s < 0 )  { s = 0.0; }
1055          zi = p.z() + s*v.z() ;
1056          if ( std::fabs(zi) <= tolODz )
1057          {
1058            xi   = p.x() + s*v.x() ;
1059            yi   = p.y() + s*v.y() ;
1060            rho2 = xi*xi + yi*yi ;
1061
1062            if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1063              || ( (rho2 >  tolORMin2) && (rho2 <  tolIRMin2)
1064                && ( v.y()*cosSPhi - v.x()*sinSPhi >  0 )
1065                && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 )     )
1066              || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1067                && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1068                && (v.x()*cosSPhi + v.y()*sinSPhi < 0) )    )
1069            {
1070              // z and r intersections good
1071              // - check intersecting with correct half-plane
1072              //
1073              if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = s; }
1074            }
1075          }
1076        }
1077      }   
1078    }
1079     
1080    // Second phi surface (Ending phi)
1081
1082    Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1083       
1084    if (Comp < 0 )  // Component in outwards normal dirn
1085    {
1086      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1087
1088      if ( Dist < halfCarTolerance )
1089      {
1090        s = Dist/Comp ;
1091
1092        if (s < snxt)
1093        {
1094          if ( s < 0 )  { s = 0; }
1095          zi = p.z() + s*v.z() ;
1096          if ( std::fabs(zi) <= tolODz )
1097          {
1098            xi   = p.x() + s*v.x() ;
1099            yi   = p.y() + s*v.y() ;
1100            rho2 = xi*xi + yi*yi ;
1101            if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1102                || ( (rho2 > tolORMin2)  && (rho2 < tolIRMin2)
1103                  && (v.x()*sinEPhi - v.y()*cosEPhi >  0)
1104                  && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1105                || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1106                  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1107                  && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1108            {
1109              // z and r intersections good
1110              // - check intersecting with correct half-plane
1111              //
1112              if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = s; }
1113            }                         //?? >=-halfCarTolerance
1114          }
1115        }
1116      }
1117    }         //  Comp < 0
1118  }           //  !fPhiFullTube
1119  if ( snxt<halfCarTolerance )  { snxt=0; }
1120  return snxt ;
1121}
1122 
1123//////////////////////////////////////////////////////////////////
1124//
1125// Calculate distance to shape from outside, along normalised vector
1126// - return kInfinity if no intersection, or intersection distance <= tolerance
1127//
1128// - Compute the intersection with the z planes
1129//        - if at valid r, phi, return
1130//
1131// -> If point is outer outer radius, compute intersection with rmax
1132//        - if at valid phi,z return
1133//
1134// -> Compute intersection with inner radius, taking largest +ve root
1135//        - if valid (in z,phi), save intersction
1136//
1137//    -> If phi segmented, compute intersections with phi half planes
1138//        - return smallest of valid phi intersections and
1139//          inner radius intersection
1140//
1141// NOTE:
1142// - Precalculations for phi trigonometry are Done `just in time'
1143// - `if valid' implies tolerant checking of intersection points
1144//   Calculate distance (<= actual) to closest surface of shape from outside
1145// - Calculate distance to z, radial planes
1146// - Only to phi planes if outside phi extent
1147// - Return 0 if point inside
1148
1149G4double G4Tubs::DistanceToIn( const G4ThreeVector& p ) const
1150{
1151  G4double safe=0.0, rho, safe1, safe2, safe3 ;
1152  G4double safePhi, cosPsi ;
1153
1154  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1155  safe1 = fRMin - rho ;
1156  safe2 = rho - fRMax ;
1157  safe3 = std::fabs(p.z()) - fDz ;
1158
1159  if ( safe1 > safe2 ) { safe = safe1; }
1160  else                 { safe = safe2; }
1161  if ( safe3 > safe )  { safe = safe3; }
1162
1163  if ( (!fPhiFullTube) && (rho) )
1164  {
1165    // Psi=angle from central phi to point
1166    //
1167    cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1168   
1169    if ( cosPsi < std::cos(fDPhi*0.5) )
1170    {
1171      // Point lies outside phi range
1172
1173      if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1174      {
1175        safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1176      }
1177      else
1178      {
1179        safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1180      }
1181      if ( safePhi > safe )  { safe = safePhi; }
1182    }
1183  }
1184  if ( safe < 0 )  { safe = 0; }
1185  return safe ;
1186}
1187
1188//////////////////////////////////////////////////////////////////////////////
1189//
1190// Calculate distance to surface of shape from `inside', allowing for tolerance
1191// - Only Calc rmax intersection if no valid rmin intersection
1192
1193G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1194                                const G4ThreeVector& v,
1195                                const G4bool calcNorm,
1196                                      G4bool *validNorm,
1197                                      G4ThreeVector *n    ) const
1198{ 
1199  ESide side=kNull , sider=kNull, sidephi=kNull ;
1200  G4double snxt, sr=kInfinity, sphi=kInfinity, pdist ;
1201  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1202
1203  static const G4double halfCarTolerance = kCarTolerance*0.5;
1204  static const G4double halfAngTolerance = kAngTolerance*0.5;
1205 
1206  // Vars for phi intersection:
1207
1208  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1209 
1210  // Z plane intersection
1211
1212  if (v.z() > 0 )
1213  {
1214    pdist = fDz - p.z() ;
1215    if ( pdist > halfCarTolerance )
1216    {
1217      snxt = pdist/v.z() ;
1218      side = kPZ ;
1219    }
1220    else
1221    {
1222      if (calcNorm)
1223      {
1224        *n         = G4ThreeVector(0,0,1) ;
1225        *validNorm = true ;
1226      }
1227      return snxt = 0 ;
1228    }
1229  }
1230  else if ( v.z() < 0 )
1231  {
1232    pdist = fDz + p.z() ;
1233
1234    if ( pdist > halfCarTolerance )
1235    {
1236      snxt = -pdist/v.z() ;
1237      side = kMZ ;
1238    }
1239    else
1240    {
1241      if (calcNorm)
1242      {
1243        *n         = G4ThreeVector(0,0,-1) ;
1244        *validNorm = true ;
1245      }
1246      return snxt = 0.0 ;
1247    }
1248  }
1249  else
1250  {
1251    snxt = kInfinity ;    // Travel perpendicular to z axis
1252    side = kNull;
1253  }
1254
1255  // Radial Intersections
1256  //
1257  // Find intersection with cylinders at rmax/rmin
1258  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1259  //
1260  // Intersects with x^2+y^2=R^2
1261  //
1262  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1263  //
1264  //            t1                t2                    t3
1265
1266  t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1267  t2   = p.x()*v.x() + p.y()*v.y() ;
1268  t3   = p.x()*p.x() + p.y()*p.y() ;
1269
1270  if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fRMax*fRMax; }
1271  else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }        // radius^2 on +-fDz
1272
1273  if ( t1 > 0 ) // Check not parallel
1274  {
1275    // Calculate sr, r exit distance
1276     
1277    if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1278    {
1279      // Delta r not negative => leaving via rmax
1280
1281      deltaR = t3 - fRMax*fRMax ;
1282
1283      // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1284      // - avoid sqrt for efficiency
1285
1286      if ( deltaR < -kRadTolerance*fRMax )
1287      {
1288        b     = t2/t1 ;
1289        c     = deltaR/t1 ;
1290        d2    = b*b-c;
1291        if( d2 >= 0 ) { sr = -b + std::sqrt(d2); }
1292        else          { sr = 0.; }
1293        sider = kRMax ;
1294      }
1295      else
1296      {
1297        // On tolerant boundary & heading outwards (or perpendicular to)
1298        // outer radial surface -> leaving immediately
1299
1300        if ( calcNorm ) 
1301        {
1302          *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1303          *validNorm = true ;
1304        }
1305        return snxt = 0 ; // Leaving by rmax immediately
1306      }
1307    }             
1308    else if ( t2 < 0. ) // i.e.  t2 < 0; Possible rmin intersection
1309    {
1310      roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1311
1312      if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1313      {
1314        deltaR = t3 - fRMin*fRMin ;
1315        b      = t2/t1 ;
1316        c      = deltaR/t1 ;
1317        d2     = b*b - c ;
1318
1319        if ( d2 >= 0 )   // Leaving via rmin
1320        {
1321          // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1322          // - avoid sqrt for efficiency
1323
1324          if (deltaR > kRadTolerance*fRMin)
1325          {
1326            sr    = -b-std::sqrt(d2) ;
1327            sider = kRMin ;
1328          }
1329          else
1330          {
1331            if ( calcNorm ) { *validNorm = false; }  // Concave side
1332            return snxt = 0.0;
1333          }
1334        }
1335        else    // No rmin intersect -> must be rmax intersect
1336        {
1337          deltaR = t3 - fRMax*fRMax ;
1338          c     = deltaR/t1 ;
1339          d2    = b*b-c;
1340          if( d2 >=0. )
1341          {
1342            sr     = -b + std::sqrt(d2) ;
1343            sider  = kRMax ;
1344          }
1345          else // Case: On the border+t2<kRadTolerance
1346               //       (v is perpendicular to the surface)
1347          {
1348            if (calcNorm)
1349            {
1350              *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1351              *validNorm = true ;
1352            }
1353            return snxt = 0.0;
1354          }
1355        }
1356      }
1357      else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1358           // No rmin intersect -> must be rmax intersect
1359      {
1360        deltaR = t3 - fRMax*fRMax ;
1361        b      = t2/t1 ;
1362        c      = deltaR/t1;
1363        d2     = b*b-c;
1364        if( d2 >= 0 )
1365        {
1366          sr     = -b + std::sqrt(d2) ;
1367          sider  = kRMax ;
1368        }
1369        else // Case: On the border+t2<kRadTolerance
1370             //       (v is perpendicular to the surface)
1371        {
1372          if (calcNorm)
1373          {
1374            *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1375            *validNorm = true ;
1376          }
1377          return snxt = 0.0;
1378        }
1379      }
1380    }
1381   
1382    // Phi Intersection
1383
1384    if ( !fPhiFullTube )
1385    {
1386      // add angle calculation with correction
1387      // of the difference in domain of atan2 and Sphi
1388      //
1389      vphi = std::atan2(v.y(),v.x()) ;
1390     
1391      if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
1392      else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1393
1394
1395      if ( p.x() || p.y() )  // Check if on z axis (rho not needed later)
1396      {
1397        // pDist -ve when inside
1398
1399        pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1400        pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1401
1402        // Comp -ve when in direction of outwards normal
1403
1404        compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1405        compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
1406       
1407        sidephi = kNull;
1408       
1409        if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1410                              && (pDistE <= halfCarTolerance) ) )
1411         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
1412                              && (pDistE >  halfCarTolerance) ) )  )
1413        {
1414          // Inside both phi *full* planes
1415         
1416          if ( compS < 0 )
1417          {
1418            sphi = pDistS/compS ;
1419           
1420            if (sphi >= -halfCarTolerance)
1421            {
1422              xi = p.x() + sphi*v.x() ;
1423              yi = p.y() + sphi*v.y() ;
1424             
1425              // Check intersecting with correct half-plane
1426              // (if not -> no intersect)
1427              //
1428              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
1429              {
1430                sidephi = kSPhi;
1431                if (((fSPhi-halfAngTolerance)<=vphi)
1432                   &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1433                {
1434                  sphi = kInfinity;
1435                }
1436              }
1437              else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1438              {
1439                sphi = kInfinity ;
1440              }
1441              else
1442              {
1443                sidephi = kSPhi ;
1444                if ( pDistS > -halfCarTolerance )
1445                {
1446                  sphi = 0.0 ; // Leave by sphi immediately
1447                }   
1448              }       
1449            }
1450            else
1451            {
1452              sphi = kInfinity ;
1453            }
1454          }
1455          else
1456          {
1457            sphi = kInfinity ;
1458          }
1459
1460          if ( compE < 0 )
1461          {
1462            sphi2 = pDistE/compE ;
1463           
1464            // Only check further if < starting phi intersection
1465            //
1466            if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1467            {
1468              xi = p.x() + sphi2*v.x() ;
1469              yi = p.y() + sphi2*v.y() ;
1470             
1471              if ((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
1472              {
1473                // Leaving via ending phi
1474                //
1475                if( !((fSPhi-halfAngTolerance <= vphi)
1476                     &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1477                {
1478                  sidephi = kEPhi ;
1479                  if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
1480                  else                                { sphi = 0.0 ;   }
1481                }
1482              } 
1483              else    // Check intersecting with correct half-plane
1484
1485              if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1486              {
1487                // Leaving via ending phi
1488                //
1489                sidephi = kEPhi ;
1490                if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1491                else                               { sphi = 0.0 ;   }
1492              }
1493            }
1494          }
1495        }
1496        else
1497        {
1498          sphi = kInfinity ;
1499        }
1500      }
1501      else
1502      {
1503        // On z axis + travel not || to z axis -> if phi of vector direction
1504        // within phi of shape, Step limited by rmax, else Step =0
1505               
1506        if ( (fSPhi - halfAngTolerance <= vphi)
1507           && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1508        {
1509          sphi = kInfinity ;
1510        }
1511        else
1512        {
1513          sidephi = kSPhi ; // arbitrary
1514          sphi    = 0.0 ;
1515        }
1516      }
1517      if (sphi < snxt)  // Order intersecttions
1518      {
1519        snxt = sphi ;
1520        side = sidephi ;
1521      }
1522    }
1523    if (sr < snxt)  // Order intersections
1524    {
1525      snxt = sr ;
1526      side = sider ;
1527    }
1528  }
1529  if (calcNorm)
1530  {
1531    switch(side)
1532    {
1533      case kRMax:
1534        // Note: returned vector not normalised
1535        // (divide by fRMax for unit vector)
1536        //
1537        xi = p.x() + snxt*v.x() ;
1538        yi = p.y() + snxt*v.y() ;
1539        *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1540        *validNorm = true ;
1541        break ;
1542
1543      case kRMin:
1544        *validNorm = false ;  // Rmin is inconvex
1545        break ;
1546
1547      case kSPhi:
1548        if ( fDPhi <= pi )
1549        {
1550          *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1551          *validNorm = true ;
1552        }
1553        else
1554        {
1555          *validNorm = false ;
1556        }
1557        break ;
1558
1559      case kEPhi:
1560        if (fDPhi <= pi)
1561        {
1562          *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1563          *validNorm = true ;
1564        }
1565        else
1566        {
1567          *validNorm = false ;
1568        }
1569        break ;
1570
1571      case kPZ:
1572        *n         = G4ThreeVector(0,0,1) ;
1573        *validNorm = true ;
1574        break ;
1575
1576      case kMZ:
1577        *n         = G4ThreeVector(0,0,-1) ;
1578        *validNorm = true ;
1579        break ;
1580
1581      default:
1582        G4cout.precision(16) ;
1583        G4cout << G4endl ;
1584        DumpInfo();
1585        G4cout << "Position:"  << G4endl << G4endl ;
1586        G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1587        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1588        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1589        G4cout << "Direction:" << G4endl << G4endl ;
1590        G4cout << "v.x() = "   << v.x() << G4endl ;
1591        G4cout << "v.y() = "   << v.y() << G4endl ;
1592        G4cout << "v.z() = "   << v.z() << G4endl << G4endl ;
1593        G4cout << "Proposed distance :" << G4endl << G4endl ;
1594        G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl ;
1595        G4Exception("G4Tubs::DistanceToOut(p,v,..)","Notification",JustWarning,
1596                    "Undefined side for valid surface normal to solid.");
1597        break ;
1598    }
1599  }
1600  if ( snxt<halfCarTolerance )  { snxt=0 ; }
1601
1602  return snxt ;
1603}
1604
1605//////////////////////////////////////////////////////////////////////////
1606//
1607// Calculate distance (<=actual) to closest surface of shape from inside
1608
1609G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
1610{
1611  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1612  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1613
1614#ifdef G4CSGDEBUG
1615  if( Inside(p) == kOutside )
1616  {
1617    G4cout.precision(16) ;
1618    G4cout << G4endl ;
1619    DumpInfo();
1620    G4cout << "Position:"  << G4endl << G4endl ;
1621    G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1622    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1623    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1624    G4Exception("G4Tubs::DistanceToOut(p)", "Notification", JustWarning, 
1625                 "Point p is outside !?");
1626  }
1627#endif
1628
1629  if ( fRMin )
1630  {
1631    safeR1 = rho   - fRMin ;
1632    safeR2 = fRMax - rho ;
1633 
1634    if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1635    else                   { safe = safeR2 ; }
1636  }
1637  else
1638  {
1639    safe = fRMax - rho ;
1640  }
1641  safeZ = fDz - std::fabs(p.z()) ;
1642
1643  if ( safeZ < safe )  { safe = safeZ ; }
1644
1645  // Check if phi divided, Calc distances closest phi plane
1646  //
1647  if ( !fPhiFullTube )
1648  {
1649    if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1650    {
1651      safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1652    }
1653    else
1654    {
1655      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1656    }
1657    if (safePhi < safe)  { safe = safePhi ; }
1658  }
1659  if ( safe < 0 )  { safe = 0 ; }
1660
1661  return safe ; 
1662}
1663
1664/////////////////////////////////////////////////////////////////////////
1665//
1666// Create a List containing the transformed vertices
1667// Ordering [0-3] -fDz cross section
1668//          [4-7] +fDz cross section such that [0] is below [4],
1669//                                             [1] below [5] etc.
1670// Note:
1671//  Caller has deletion resposibility
1672//  Potential improvement: For last slice, use actual ending angle
1673//                         to avoid rounding error problems.
1674
1675G4ThreeVectorList*
1676G4Tubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1677{
1678  G4ThreeVectorList* vertices ;
1679  G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
1680  G4double meshAngle, meshRMax, crossAngle,
1681           cosCrossAngle, sinCrossAngle, sAngle;
1682  G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1683  G4int crossSection, noCrossSections;
1684
1685  // Compute no of cross-sections necessary to mesh tube
1686  //
1687  noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
1688
1689  if ( noCrossSections < kMinMeshSections )
1690  {
1691    noCrossSections = kMinMeshSections ;
1692  }
1693  else if (noCrossSections>kMaxMeshSections)
1694  {
1695    noCrossSections = kMaxMeshSections ;
1696  }
1697  // noCrossSections = 4 ;
1698
1699  meshAngle = fDPhi/(noCrossSections - 1) ;
1700  // meshAngle = fDPhi/(noCrossSections) ;
1701
1702  meshRMax  = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
1703  meshRMin = fRMin - 100*kCarTolerance ; 
1704 
1705  // If complete in phi, set start angle such that mesh will be at fRMax
1706  // on the x axis. Will give better extent calculations when not rotated.
1707
1708  if (fPhiFullTube && (fSPhi == 0) )  { sAngle = -meshAngle*0.5 ; }
1709  else                                { sAngle =  fSPhi ; }
1710   
1711  vertices = new G4ThreeVectorList();
1712  vertices->reserve(noCrossSections*4);
1713   
1714  if ( vertices )
1715  {
1716    for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1717    {
1718      // Compute coordinates of cross section at section crossSection
1719
1720      crossAngle    = sAngle + crossSection*meshAngle ;
1721      cosCrossAngle = std::cos(crossAngle) ;
1722      sinCrossAngle = std::sin(crossAngle) ;
1723
1724      rMaxX = meshRMax*cosCrossAngle ;
1725      rMaxY = meshRMax*sinCrossAngle ;
1726
1727      if(meshRMin <= 0.0)
1728      {
1729        rMinX = 0.0 ;
1730        rMinY = 0.0 ;
1731      }
1732      else
1733      {
1734        rMinX = meshRMin*cosCrossAngle ;
1735        rMinY = meshRMin*sinCrossAngle ;
1736      }
1737      vertex0 = G4ThreeVector(rMinX,rMinY,-fDz) ;
1738      vertex1 = G4ThreeVector(rMaxX,rMaxY,-fDz) ;
1739      vertex2 = G4ThreeVector(rMaxX,rMaxY,+fDz) ;
1740      vertex3 = G4ThreeVector(rMinX,rMinY,+fDz) ;
1741
1742      vertices->push_back(pTransform.TransformPoint(vertex0)) ;
1743      vertices->push_back(pTransform.TransformPoint(vertex1)) ;
1744      vertices->push_back(pTransform.TransformPoint(vertex2)) ;
1745      vertices->push_back(pTransform.TransformPoint(vertex3)) ;
1746    }
1747  }
1748  else
1749  {
1750    DumpInfo();
1751    G4Exception("G4Tubs::CreateRotatedVertices()",
1752                "FatalError", FatalException,
1753                "Error in allocation of vertices. Out of memory !");
1754  }
1755  return vertices ;
1756}
1757
1758//////////////////////////////////////////////////////////////////////////
1759//
1760// Stream object contents to an output stream
1761
1762G4GeometryType G4Tubs::GetEntityType() const
1763{
1764  return G4String("G4Tubs");
1765}
1766
1767//////////////////////////////////////////////////////////////////////////
1768//
1769// Stream object contents to an output stream
1770
1771std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1772{
1773  os << "-----------------------------------------------------------\n"
1774     << "    *** Dump for solid - " << GetName() << " ***\n"
1775     << "    ===================================================\n"
1776     << " Solid type: G4Tubs\n"
1777     << " Parameters: \n"
1778     << "    inner radius : " << fRMin/mm << " mm \n"
1779     << "    outer radius : " << fRMax/mm << " mm \n"
1780     << "    half length Z: " << fDz/mm << " mm \n"
1781     << "    starting phi : " << fSPhi/degree << " degrees \n"
1782     << "    delta phi    : " << fDPhi/degree << " degrees \n"
1783     << "-----------------------------------------------------------\n";
1784
1785  return os;
1786}
1787
1788/////////////////////////////////////////////////////////////////////////
1789//
1790// GetPointOnSurface
1791
1792G4ThreeVector G4Tubs::GetPointOnSurface() const
1793{
1794  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1795           aOne, aTwo, aThr, aFou;
1796  G4double rRand;
1797
1798  aOne = 2.*fDz*fDPhi*fRMax;
1799  aTwo = 2.*fDz*fDPhi*fRMin;
1800  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1801  aFou = 2.*fDz*(fRMax-fRMin);
1802
1803  phi    = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1804  cosphi = std::cos(phi);
1805  sinphi = std::sin(phi);
1806
1807  rRand  = RandFlat::shoot(fRMin,fRMax);
1808 
1809  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1810 
1811  chose  = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1812
1813  if( (chose >=0) && (chose < aOne) )
1814  {
1815    xRand = fRMax*cosphi;
1816    yRand = fRMax*sinphi;
1817    zRand = RandFlat::shoot(-1.*fDz,fDz);
1818    return G4ThreeVector  (xRand, yRand, zRand);
1819  }
1820  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1821  {
1822    xRand = fRMin*cosphi;
1823    yRand = fRMin*sinphi;
1824    zRand = RandFlat::shoot(-1.*fDz,fDz);
1825    return G4ThreeVector  (xRand, yRand, zRand);
1826  }
1827  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1828  {
1829    xRand = rRand*cosphi;
1830    yRand = rRand*sinphi;
1831    zRand = fDz;
1832    return G4ThreeVector  (xRand, yRand, zRand);
1833  }
1834  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1835  {
1836    xRand = rRand*cosphi;
1837    yRand = rRand*sinphi;
1838    zRand = -1.*fDz;
1839    return G4ThreeVector  (xRand, yRand, zRand);
1840  }
1841  else if( (chose >= aOne + aTwo + 2.*aThr)
1842        && (chose < aOne + aTwo + 2.*aThr + aFou) )
1843  {
1844    xRand = rRand*std::cos(fSPhi);
1845    yRand = rRand*std::sin(fSPhi);
1846    zRand = RandFlat::shoot(-1.*fDz,fDz);
1847    return G4ThreeVector  (xRand, yRand, zRand);
1848  }
1849  else
1850  {
1851    xRand = rRand*std::cos(fSPhi+fDPhi);
1852    yRand = rRand*std::sin(fSPhi+fDPhi);
1853    zRand = RandFlat::shoot(-1.*fDz,fDz);
1854    return G4ThreeVector  (xRand, yRand, zRand);
1855  }
1856}
1857
1858///////////////////////////////////////////////////////////////////////////
1859//
1860// Methods for visualisation
1861
1862void G4Tubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
1863{
1864  scene.AddSolid (*this) ;
1865}
1866
1867G4Polyhedron* G4Tubs::CreatePolyhedron () const 
1868{
1869  return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1870}
1871
1872G4NURBS* G4Tubs::CreateNURBS () const 
1873{
1874  G4NURBS* pNURBS ;
1875  if (fRMin != 0) 
1876  {
1877    if (fPhiFullTube) 
1878    {
1879      pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
1880    }
1881    else 
1882    {
1883      pNURBS = new G4NURBStubesector (fRMin,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1884    }
1885  }
1886  else 
1887  {
1888    if (fPhiFullTube) 
1889    {
1890      pNURBS = new G4NURBScylinder (fRMax,fDz) ;
1891    }
1892    else 
1893    {
1894      const G4double epsilon = 1.e-4 ; // Cylinder sector not yet available!
1895      pNURBS = new G4NURBStubesector (epsilon,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1896    }
1897  }
1898  return pNURBS ;
1899}
Note: See TracBrowser for help on using the repository browser.