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

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

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