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