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

Last change on this file since 836 was 831, checked in by garnier, 16 years ago

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