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