Geant4.10
 All Data Structures 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 76732 2013-11-14 14:36: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 #if !defined(G4GEOM_USE_UCONS)
48 
49 #include "G4VoxelLimits.hh"
50 #include "G4AffineTransform.hh"
51 #include "G4GeometryTolerance.hh"
52 
53 #include "G4VPVParameterisation.hh"
54 
55 #include "meshdefs.hh"
56 
57 #include "Randomize.hh"
58 
59 #include "G4VGraphicsScene.hh"
60 #include "G4Polyhedron.hh"
61 
62 using namespace CLHEP;
63 
64 ////////////////////////////////////////////////////////////////////////
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 
74 //////////////////////////////////////////////////////////////////////////
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  halfCarTolerance=kCarTolerance*0.5;
91  halfRadTolerance=kRadTolerance*0.5;
92  halfAngTolerance=kAngTolerance*0.5;
93 
94  // Check z-len
95  //
96  if ( pDz < 0 )
97  {
98  std::ostringstream message;
99  message << "Invalid Z half-length for Solid: " << GetName() << G4endl
100  << " hZ = " << pDz;
101  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
102  FatalException, message);
103  }
104 
105  // Check radii
106  //
107  if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
108  {
109  std::ostringstream message;
110  message << "Invalid values of radii for Solid: " << GetName() << G4endl
111  << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
112  << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
113  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
114  FatalException, message) ;
115  }
116  if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
117  if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
118 
119  // Check angles
120  //
121  CheckPhiAngles(pSPhi, pDPhi);
122 }
123 
124 ///////////////////////////////////////////////////////////////////////
125 //
126 // Fake default constructor - sets only member data and allocates memory
127 // for usage restricted to object persistency.
128 //
129 G4Cons::G4Cons( __void__& a )
130  : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
131  fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
132  fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
133  cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
134  fPhiFullCone(false), halfCarTolerance(0.), halfRadTolerance(0.),
135  halfAngTolerance(0.)
136 
137 {
138 }
139 
140 ///////////////////////////////////////////////////////////////////////
141 //
142 // Destructor
143 
144 G4Cons::~G4Cons()
145 {
146 }
147 
148 //////////////////////////////////////////////////////////////////////////
149 //
150 // Copy constructor
151 
152 G4Cons::G4Cons(const G4Cons& rhs)
153  : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
154  kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
155  fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
156  fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
157  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
158  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
159  cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
160  halfCarTolerance(rhs.halfCarTolerance),
161  halfRadTolerance(rhs.halfRadTolerance),
162  halfAngTolerance(rhs.halfAngTolerance)
163 {
164 }
165 
166 //////////////////////////////////////////////////////////////////////////
167 //
168 // Assignment operator
169 
170 G4Cons& G4Cons::operator = (const G4Cons& rhs)
171 {
172  // Check assignment to self
173  //
174  if (this == &rhs) { return *this; }
175 
176  // Copy base class data
177  //
178  G4CSGSolid::operator=(rhs);
179 
180  // Copy data
181  //
182  kRadTolerance = rhs.kRadTolerance;
183  kAngTolerance = rhs.kAngTolerance;
184  fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
185  fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
186  fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
187  sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
188  cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
189  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
190  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
191  fPhiFullCone = rhs.fPhiFullCone;
192  halfCarTolerance = rhs.halfCarTolerance;
193  halfRadTolerance = rhs.halfRadTolerance;
194  halfAngTolerance = rhs.halfAngTolerance;
195 
196  return *this;
197 }
198 
199 /////////////////////////////////////////////////////////////////////
200 //
201 // Return whether point inside/outside/on surface
202 
203 EInside G4Cons::Inside(const G4ThreeVector& p) const
204 {
205  G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
206  EInside in;
207 
208  if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
209  else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
210  else { in = kInside; }
211 
212  r2 = p.x()*p.x() + p.y()*p.y() ;
213  rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
214  rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
215 
216  // rh2 = rh*rh;
217 
218  tolRMin = rl - halfRadTolerance;
219  if ( tolRMin < 0 ) { tolRMin = 0; }
220  tolRMax = rh + halfRadTolerance;
221 
222  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
223 
224  if (rl) { tolRMin = rl + halfRadTolerance; }
225  else { tolRMin = 0.0; }
226  tolRMax = rh - halfRadTolerance;
227 
228  if (in == kInside) // else it's kSurface already
229  {
230  if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
231  }
232  if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
233  {
234  pPhi = std::atan2(p.y(),p.x()) ;
235 
236  if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
237  else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
238 
239  if ( (pPhi < fSPhi - halfAngTolerance) ||
240  (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
241 
242  else if (in == kInside) // else it's kSurface anyway already
243  {
244  if ( (pPhi < fSPhi + halfAngTolerance) ||
245  (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
246  }
247  }
248  else if ( !fPhiFullCone ) { in = kSurface; }
249 
250  return in ;
251 }
252 
253 /////////////////////////////////////////////////////////////////////////
254 //
255 // Dispatch to parameterisation for replication mechanism dimension
256 // computation & modification.
257 
258 void G4Cons::ComputeDimensions( G4VPVParameterisation* p,
259  const G4int n,
260  const G4VPhysicalVolume* pRep )
261 {
262  p->ComputeDimensions(*this,n,pRep) ;
263 }
264 
265 
266 ///////////////////////////////////////////////////////////////////////////
267 //
268 // Calculate extent under transform and specified limit
269 
270 G4bool G4Cons::CalculateExtent( const EAxis pAxis,
271  const G4VoxelLimits& pVoxelLimit,
272  const G4AffineTransform& pTransform,
273  G4double& pMin,
274  G4double& pMax ) const
275 {
276  if ( !pTransform.IsRotated() && (fDPhi == twopi)
277  && (fRmin1 == 0) && (fRmin2 == 0) )
278  {
279  // Special case handling for unrotated solid cones
280  // Compute z/x/y mins and maxs for bounding box respecting limits,
281  // with early returns if outside limits. Then switch() on pAxis,
282  // and compute exact x and y limit for x/y case
283 
284  G4double xoffset, xMin, xMax ;
285  G4double yoffset, yMin, yMax ;
286  G4double zoffset, zMin, zMax ;
287 
288  G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
289  G4double xoff1, xoff2, yoff1, yoff2 ;
290 
291  zoffset = pTransform.NetTranslation().z();
292  zMin = zoffset - fDz ;
293  zMax = zoffset + fDz ;
294 
295  if (pVoxelLimit.IsZLimited())
296  {
297  if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
298  (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
299  {
300  return false ;
301  }
302  else
303  {
304  if ( zMin < pVoxelLimit.GetMinZExtent() )
305  {
306  zMin = pVoxelLimit.GetMinZExtent() ;
307  }
308  if ( zMax > pVoxelLimit.GetMaxZExtent() )
309  {
310  zMax = pVoxelLimit.GetMaxZExtent() ;
311  }
312  }
313  }
314  xoffset = pTransform.NetTranslation().x() ;
315  RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
316  xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
317  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
318  xMin = 2*xoffset-xMax ;
319 
320  if (pVoxelLimit.IsXLimited())
321  {
322  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
323  (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
324  {
325  return false ;
326  }
327  else
328  {
329  if ( xMin < pVoxelLimit.GetMinXExtent() )
330  {
331  xMin = pVoxelLimit.GetMinXExtent() ;
332  }
333  if ( xMax > pVoxelLimit.GetMaxXExtent() )
334  {
335  xMax=pVoxelLimit.GetMaxXExtent() ;
336  }
337  }
338  }
339  yoffset = pTransform.NetTranslation().y() ;
340  yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
341  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
342  yMin = 2*yoffset-yMax ;
343  RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings
344 
345  if (pVoxelLimit.IsYLimited())
346  {
347  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
348  (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
349  {
350  return false ;
351  }
352  else
353  {
354  if ( yMin < pVoxelLimit.GetMinYExtent() )
355  {
356  yMin = pVoxelLimit.GetMinYExtent() ;
357  }
358  if ( yMax > pVoxelLimit.GetMaxYExtent() )
359  {
360  yMax = pVoxelLimit.GetMaxYExtent() ;
361  }
362  }
363  }
364  switch (pAxis) // Known to cut cones
365  {
366  case kXAxis:
367  yoff1 = yoffset - yMin ;
368  yoff2 = yMax - yoffset ;
369 
370  if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
371  { // => no change
372  pMin = xMin ;
373  pMax = xMax ;
374  }
375  else
376  {
377  // Y limits don't cross max/min x => compute max delta x,
378  // hence new mins/maxs
379  delta=RMax*RMax-yoff1*yoff1;
380  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
381  delta=RMax*RMax-yoff2*yoff2;
382  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
383  maxDiff = (diff1>diff2) ? diff1:diff2 ;
384  newMin = xoffset - maxDiff ;
385  newMax = xoffset + maxDiff ;
386  pMin = ( newMin < xMin ) ? xMin : newMin ;
387  pMax = ( newMax > xMax) ? xMax : newMax ;
388  }
389  break ;
390 
391  case kYAxis:
392  xoff1 = xoffset - xMin ;
393  xoff2 = xMax - xoffset ;
394 
395  if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
396  { // => no change
397  pMin = yMin ;
398  pMax = yMax ;
399  }
400  else
401  {
402  // X limits don't cross max/min y => compute max delta y,
403  // hence new mins/maxs
404  delta=RMax*RMax-xoff1*xoff1;
405  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
406  delta=RMax*RMax-xoff2*xoff2;
407  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
408  maxDiff = (diff1 > diff2) ? diff1:diff2 ;
409  newMin = yoffset - maxDiff ;
410  newMax = yoffset + maxDiff ;
411  pMin = (newMin < yMin) ? yMin : newMin ;
412  pMax = (newMax > yMax) ? yMax : newMax ;
413  }
414  break ;
415 
416  case kZAxis:
417  pMin = zMin ;
418  pMax = zMax ;
419  break ;
420 
421  default:
422  break ;
423  }
424  pMin -= kCarTolerance ;
425  pMax += kCarTolerance ;
426 
427  return true ;
428  }
429  else // Calculate rotated vertex coordinates
430  {
431  G4int i, noEntries, noBetweenSections4 ;
432  G4bool existsAfterClip = false ;
433  G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
434 
435  pMin = +kInfinity ;
436  pMax = -kInfinity ;
437 
438  noEntries = vertices->size() ;
439  noBetweenSections4 = noEntries-4 ;
440 
441  for ( i = 0 ; i < noEntries ; i += 4 )
442  {
443  ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
444  }
445  for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
446  {
447  ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
448  }
449  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
450  {
451  existsAfterClip = true ;
452 
453  // Add 2*tolerance to avoid precision troubles
454 
455  pMin -= kCarTolerance ;
456  pMax += kCarTolerance ;
457  }
458  else
459  {
460  // Check for case where completely enveloping clipping volume
461  // If point inside then we are confident that the solid completely
462  // envelopes the clipping volume. Hence set min/max extents according
463  // to clipping volume extents along the specified axis.
464 
465  G4ThreeVector clipCentre(
466  (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
467  (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
468  (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ;
469 
470  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
471  {
472  existsAfterClip = true ;
473  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
474  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
475  }
476  }
477  delete vertices ;
478  return existsAfterClip ;
479  }
480 }
481 
482 ////////////////////////////////////////////////////////////////////////
483 //
484 // Return unit normal of surface closest to p
485 // - note if point on z axis, ignore phi divided sides
486 // - unsafe if point close to z axis a rmin=0 - no explicit checks
487 
488 G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const
489 {
490  G4int noSurfaces = 0;
491  G4double rho, pPhi;
492  G4double distZ, distRMin, distRMax;
493  G4double distSPhi = kInfinity, distEPhi = kInfinity;
494  G4double tanRMin, secRMin, pRMin, widRMin;
495  G4double tanRMax, secRMax, pRMax, widRMax;
496 
497  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
498  G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
499 
500  distZ = std::fabs(std::fabs(p.z()) - fDz);
501  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
502 
503  tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
504  secRMin = std::sqrt(1 + tanRMin*tanRMin);
505  pRMin = rho - p.z()*tanRMin;
506  widRMin = fRmin2 - fDz*tanRMin;
507  distRMin = std::fabs(pRMin - widRMin)/secRMin;
508 
509  tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
510  secRMax = std::sqrt(1+tanRMax*tanRMax);
511  pRMax = rho - p.z()*tanRMax;
512  widRMax = fRmax2 - fDz*tanRMax;
513  distRMax = std::fabs(pRMax - widRMax)/secRMax;
514 
515  if (!fPhiFullCone) // Protected against (0,0,z)
516  {
517  if ( rho )
518  {
519  pPhi = std::atan2(p.y(),p.x());
520 
521  if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
522  else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
523 
524  distSPhi = std::fabs( pPhi - fSPhi );
525  distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
526  }
527  else if( !(fRmin1) || !(fRmin2) )
528  {
529  distSPhi = 0.;
530  distEPhi = 0.;
531  }
532  nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
533  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
534  }
535  if ( rho > halfCarTolerance )
536  {
537  nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
538  if (fRmin1 || fRmin2)
539  {
540  nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
541  }
542  }
543 
544  if( distRMax <= halfCarTolerance )
545  {
546  noSurfaces ++;
547  sumnorm += nR;
548  }
549  if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
550  {
551  noSurfaces ++;
552  sumnorm += nr;
553  }
554  if( !fPhiFullCone )
555  {
556  if (distSPhi <= halfAngTolerance)
557  {
558  noSurfaces ++;
559  sumnorm += nPs;
560  }
561  if (distEPhi <= halfAngTolerance)
562  {
563  noSurfaces ++;
564  sumnorm += nPe;
565  }
566  }
567  if (distZ <= halfCarTolerance)
568  {
569  noSurfaces ++;
570  if ( p.z() >= 0.) { sumnorm += nZ; }
571  else { sumnorm -= nZ; }
572  }
573  if ( noSurfaces == 0 )
574  {
575 #ifdef G4CSGDEBUG
576  G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
577  JustWarning, "Point p is not on surface !?" );
578 #endif
579  norm = ApproxSurfaceNormal(p);
580  }
581  else if ( noSurfaces == 1 ) { norm = sumnorm; }
582  else { norm = sumnorm.unit(); }
583 
584  return norm ;
585 }
586 
587 ////////////////////////////////////////////////////////////////////////////
588 //
589 // Algorithm for SurfaceNormal() following the original specification
590 // for points not on the surface
591 
592 G4ThreeVector G4Cons::ApproxSurfaceNormal( const G4ThreeVector& p ) const
593 {
594  ENorm side ;
595  G4ThreeVector norm ;
596  G4double rho, phi ;
597  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
598  G4double tanRMin, secRMin, pRMin, widRMin ;
599  G4double tanRMax, secRMax, pRMax, widRMax ;
600 
601  distZ = std::fabs(std::fabs(p.z()) - fDz) ;
602  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
603 
604  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
605  secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
606  pRMin = rho - p.z()*tanRMin ;
607  widRMin = fRmin2 - fDz*tanRMin ;
608  distRMin = std::fabs(pRMin - widRMin)/secRMin ;
609 
610  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
611  secRMax = std::sqrt(1+tanRMax*tanRMax) ;
612  pRMax = rho - p.z()*tanRMax ;
613  widRMax = fRmax2 - fDz*tanRMax ;
614  distRMax = std::fabs(pRMax - widRMax)/secRMax ;
615 
616  if (distRMin < distRMax) // First minimum
617  {
618  if (distZ < distRMin)
619  {
620  distMin = distZ ;
621  side = kNZ ;
622  }
623  else
624  {
625  distMin = distRMin ;
626  side = kNRMin ;
627  }
628  }
629  else
630  {
631  if (distZ < distRMax)
632  {
633  distMin = distZ ;
634  side = kNZ ;
635  }
636  else
637  {
638  distMin = distRMax ;
639  side = kNRMax ;
640  }
641  }
642  if ( !fPhiFullCone && rho ) // Protected against (0,0,z)
643  {
644  phi = std::atan2(p.y(),p.x()) ;
645 
646  if (phi < 0) { phi += twopi; }
647 
648  if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
649  else { distSPhi = std::fabs(phi - fSPhi)*rho; }
650 
651  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
652 
653  // Find new minimum
654 
655  if (distSPhi < distEPhi)
656  {
657  if (distSPhi < distMin) { side = kNSPhi; }
658  }
659  else
660  {
661  if (distEPhi < distMin) { side = kNEPhi; }
662  }
663  }
664  switch (side)
665  {
666  case kNRMin: // Inner radius
667  rho *= secRMin ;
668  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
669  break ;
670  case kNRMax: // Outer radius
671  rho *= secRMax ;
672  norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
673  break ;
674  case kNZ: // +/- dz
675  if (p.z() > 0) { norm = G4ThreeVector(0,0,1); }
676  else { norm = G4ThreeVector(0,0,-1); }
677  break ;
678  case kNSPhi:
679  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
680  break ;
681  case kNEPhi:
682  norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
683  break ;
684  default: // Should never reach this case...
685  DumpInfo();
686  G4Exception("G4Cons::ApproxSurfaceNormal()",
687  "GeomSolids1002", JustWarning,
688  "Undefined side for valid surface normal to solid.");
689  break ;
690  }
691  return norm ;
692 }
693 
694 ////////////////////////////////////////////////////////////////////////
695 //
696 // Calculate distance to shape from outside, along normalised vector
697 // - return kInfinity if no intersection, or intersection distance <= tolerance
698 //
699 // - Compute the intersection with the z planes
700 // - if at valid r, phi, return
701 //
702 // -> If point is outside cone, compute intersection with rmax1*0.5
703 // - if at valid phi,z return
704 // - if inside outer cone, handle case when on tolerant outer cone
705 // boundary and heading inwards(->0 to in)
706 //
707 // -> Compute intersection with inner cone, taking largest +ve root
708 // - if valid (in z,phi), save intersction
709 //
710 // -> If phi segmented, compute intersections with phi half planes
711 // - return smallest of valid phi intersections and
712 // inner radius intersection
713 //
714 // NOTE:
715 // - `if valid' implies tolerant checking of intersection points
716 // - z, phi intersection from Tubs
717 
718 G4double G4Cons::DistanceToIn( const G4ThreeVector& p,
719  const G4ThreeVector& v ) const
720 {
721  G4double snxt = kInfinity ; // snxt = default return value
722  const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.
723 
724  G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones
725  G4double tanRMin,secRMin,rMinAv,rMinOAv ;
726  G4double rout,rin ;
727 
728  G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared
729  G4double tolORMax2,tolIRMax,tolIRMax2 ;
730  G4double tolODz,tolIDz ;
731 
732  G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
733 
734  G4double t1,t2,t3,b,c,d ; // Quadratic solver variables
735  G4double nt1,nt2,nt3 ;
736  G4double Comp ;
737 
738  G4ThreeVector Normal;
739 
740  // Cone Precalcs
741 
742  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
743  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
744  rMinAv = (fRmin1 + fRmin2)*0.5 ;
745 
746  if (rMinAv > halfRadTolerance)
747  {
748  rMinOAv = rMinAv - halfRadTolerance ;
749  }
750  else
751  {
752  rMinOAv = 0.0 ;
753  }
754  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
755  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
756  rMaxAv = (fRmax1 + fRmax2)*0.5 ;
757  rMaxOAv = rMaxAv + halfRadTolerance ;
758 
759  // Intersection with z-surfaces
760 
761  tolIDz = fDz - halfCarTolerance ;
762  tolODz = fDz + halfCarTolerance ;
763 
764  if (std::fabs(p.z()) >= tolIDz)
765  {
766  if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa
767  {
768  sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
769 
770  if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
771 
772  xi = p.x() + sd*v.x() ; // Intersection coords
773  yi = p.y() + sd*v.y() ;
774  rhoi2 = xi*xi + yi*yi ;
775 
776  // Check validity of intersection
777  // Calculate (outer) tolerant radi^2 at intersecion
778 
779  if (v.z() > 0)
780  {
781  tolORMin = fRmin1 - halfRadTolerance*secRMin ;
782  tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
783  tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
784  tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
785  (fRmax1 + halfRadTolerance*secRMax) ;
786  }
787  else
788  {
789  tolORMin = fRmin2 - halfRadTolerance*secRMin ;
790  tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
791  tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
792  tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
793  (fRmax2 + halfRadTolerance*secRMax) ;
794  }
795  if ( tolORMin > 0 )
796  {
797  tolORMin2 = tolORMin*tolORMin ;
798  tolIRMin2 = tolIRMin*tolIRMin ;
799  }
800  else
801  {
802  tolORMin2 = 0.0 ;
803  tolIRMin2 = 0.0 ;
804  }
805  if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
806  else { tolIRMax2 = 0.0; }
807 
808  if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
809  {
810  if ( !fPhiFullCone && rhoi2 )
811  {
812  // Psi = angle made with central (average) phi of shape
813 
814  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
815 
816  if (cosPsi >= cosHDPhiIT) { return sd; }
817  }
818  else
819  {
820  return sd;
821  }
822  }
823  }
824  else // On/outside extent, and heading away -> cannot intersect
825  {
826  return snxt ;
827  }
828  }
829 
830 // ----> Can not intersect z surfaces
831 
832 
833 // Intersection with outer cone (possible return) and
834 // inner cone (must also check phi)
835 //
836 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
837 //
838 // Intersects with x^2+y^2=(a*z+b)^2
839 //
840 // where a=tanRMax or tanRMin
841 // b=rMaxAv or rMinAv
842 //
843 // (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 ;
844 // t1 t2 t3
845 //
846 // \--------u-------/ \-----------v----------/ \---------w--------/
847 //
848 
849  t1 = 1.0 - v.z()*v.z() ;
850  t2 = p.x()*v.x() + p.y()*v.y() ;
851  t3 = p.x()*p.x() + p.y()*p.y() ;
852  rin = tanRMin*p.z() + rMinAv ;
853  rout = tanRMax*p.z() + rMaxAv ;
854 
855  // Outer Cone Intersection
856  // Must be outside/on outer cone for valid intersection
857 
858  nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
859  nt2 = t2 - tanRMax*v.z()*rout ;
860  nt3 = t3 - rout*rout ;
861 
862  if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
863  {
864  b = nt2/nt1;
865  c = nt3/nt1;
866  d = b*b-c ;
867  if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
868  || (rout < 0) )
869  {
870  // If outside real cone (should be rho-rout>kRadTolerance*0.5
871  // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
872 
873  if (d >= 0)
874  {
875 
876  if ((rout < 0) && (nt3 <= 0))
877  {
878  // Inside `shadow cone' with -ve radius
879  // -> 2nd root could be on real cone
880 
881  if (b>0) { sd = c/(-b-std::sqrt(d)); }
882  else { sd = -b + std::sqrt(d); }
883  }
884  else
885  {
886  if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
887  {
888  sd=c/(-b+std::sqrt(d));
889  }
890  else
891  {
892  if ( c <= 0 ) // second >=0
893  {
894  sd = -b + std::sqrt(d) ;
895  if((sd<0) & (sd>-halfRadTolerance)) sd=0;
896  }
897  else // both negative, travel away
898  {
899  return kInfinity ;
900  }
901  }
902  }
903  if ( sd > 0 ) // If 'forwards'. Check z intersection
904  {
905  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
906  { // 64 bits systems. Split long distances and recompute
907  G4double fTerm = sd-std::fmod(sd,dRmax);
908  sd = fTerm + DistanceToIn(p+fTerm*v,v);
909  }
910  zi = p.z() + sd*v.z() ;
911 
912  if (std::fabs(zi) <= tolODz)
913  {
914  // Z ok. Check phi intersection if reqd
915 
916  if ( fPhiFullCone ) { return sd; }
917  else
918  {
919  xi = p.x() + sd*v.x() ;
920  yi = p.y() + sd*v.y() ;
921  ri = rMaxAv + zi*tanRMax ;
922  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
923 
924  if ( cosPsi >= cosHDPhiIT ) { return sd; }
925  }
926  }
927  } // end if (sd>0)
928  }
929  }
930  else
931  {
932  // Inside outer cone
933  // check not inside, and heading through G4Cons (-> 0 to in)
934 
935  if ( ( t3 > (rin + halfRadTolerance*secRMin)*
936  (rin + halfRadTolerance*secRMin) )
937  && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
938  {
939  // Inside cones, delta r -ve, inside z extent
940  // Point is on the Surface => check Direction using Normal.dot(v)
941 
942  xi = p.x() ;
943  yi = p.y() ;
944  risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
945  Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
946  if ( !fPhiFullCone )
947  {
948  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
949  if ( cosPsi >= cosHDPhiIT )
950  {
951  if ( Normal.dot(v) <= 0 ) { return 0.0; }
952  }
953  }
954  else
955  {
956  if ( Normal.dot(v) <= 0 ) { return 0.0; }
957  }
958  }
959  }
960  }
961  else // Single root case
962  {
963  if ( std::fabs(nt2) > kRadTolerance )
964  {
965  sd = -0.5*nt3/nt2 ;
966 
967  if ( sd < 0 ) { return kInfinity; } // travel away
968  else // sd >= 0, If 'forwards'. Check z intersection
969  {
970  zi = p.z() + sd*v.z() ;
971 
972  if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
973  {
974  // Z ok. Check phi intersection if reqd
975 
976  if ( fPhiFullCone ) { return sd; }
977  else
978  {
979  xi = p.x() + sd*v.x() ;
980  yi = p.y() + sd*v.y() ;
981  ri = rMaxAv + zi*tanRMax ;
982  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
983 
984  if (cosPsi >= cosHDPhiIT) { return sd; }
985  }
986  }
987  }
988  }
989  else // travel || cone surface from its origin
990  {
991  sd = kInfinity ;
992  }
993  }
994 
995  // Inner Cone Intersection
996  // o Space is divided into 3 areas:
997  // 1) Radius greater than real inner cone & imaginary cone & outside
998  // tolerance
999  // 2) Radius less than inner or imaginary cone & outside tolarance
1000  // 3) Within tolerance of real or imaginary cones
1001  // - Extra checks needed for 3's intersections
1002  // => lots of duplicated code
1003 
1004  if (rMinAv)
1005  {
1006  nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1007  nt2 = t2 - tanRMin*v.z()*rin ;
1008  nt3 = t3 - rin*rin ;
1009 
1010  if ( nt1 )
1011  {
1012  if ( nt3 > rin*kRadTolerance*secRMin )
1013  {
1014  // At radius greater than real & imaginary cones
1015  // -> 2nd root, with zi check
1016 
1017  b = nt2/nt1 ;
1018  c = nt3/nt1 ;
1019  d = b*b-c ;
1020  if (d >= 0) // > 0
1021  {
1022  if(b>0){sd = c/( -b-std::sqrt(d));}
1023  else {sd = -b + std::sqrt(d) ;}
1024 
1025  if ( sd >= 0 ) // > 0
1026  {
1027  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1028  { // 64 bits systems. Split long distance and recompute
1029  G4double fTerm = sd-std::fmod(sd,dRmax);
1030  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1031  }
1032  zi = p.z() + sd*v.z() ;
1033 
1034  if ( std::fabs(zi) <= tolODz )
1035  {
1036  if ( !fPhiFullCone )
1037  {
1038  xi = p.x() + sd*v.x() ;
1039  yi = p.y() + sd*v.y() ;
1040  ri = rMinAv + zi*tanRMin ;
1041  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1042 
1043  if (cosPsi >= cosHDPhiIT)
1044  {
1045  if ( sd > halfRadTolerance ) { snxt=sd; }
1046  else
1047  {
1048  // Calculate a normal vector in order to check Direction
1049 
1050  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1051  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1052  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1053  }
1054  }
1055  }
1056  else
1057  {
1058  if ( sd > halfRadTolerance ) { return sd; }
1059  else
1060  {
1061  // Calculate a normal vector in order to check Direction
1062 
1063  xi = p.x() + sd*v.x() ;
1064  yi = p.y() + sd*v.y() ;
1065  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1066  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1067  if ( Normal.dot(v) <= 0 ) { return sd; }
1068  }
1069  }
1070  }
1071  }
1072  }
1073  }
1074  else if ( nt3 < -rin*kRadTolerance*secRMin )
1075  {
1076  // Within radius of inner cone (real or imaginary)
1077  // -> Try 2nd root, with checking intersection is with real cone
1078  // -> If check fails, try 1st root, also checking intersection is
1079  // on real cone
1080 
1081  b = nt2/nt1 ;
1082  c = nt3/nt1 ;
1083  d = b*b - c ;
1084 
1085  if ( d >= 0 ) // > 0
1086  {
1087  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1088  else { sd = -b + std::sqrt(d); }
1089  zi = p.z() + sd*v.z() ;
1090  ri = rMinAv + zi*tanRMin ;
1091 
1092  if ( ri > 0 )
1093  {
1094  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1095  {
1096  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1097  { // seen on 64 bits systems. Split and recompute
1098  G4double fTerm = sd-std::fmod(sd,dRmax);
1099  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1100  }
1101  if ( !fPhiFullCone )
1102  {
1103  xi = p.x() + sd*v.x() ;
1104  yi = p.y() + sd*v.y() ;
1105  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1106 
1107  if (cosPsi >= cosHDPhiOT)
1108  {
1109  if ( sd > halfRadTolerance ) { snxt=sd; }
1110  else
1111  {
1112  // Calculate a normal vector in order to check Direction
1113 
1114  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1115  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1116  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1117  }
1118  }
1119  }
1120  else
1121  {
1122  if( sd > halfRadTolerance ) { return sd; }
1123  else
1124  {
1125  // Calculate a normal vector in order to check Direction
1126 
1127  xi = p.x() + sd*v.x() ;
1128  yi = p.y() + sd*v.y() ;
1129  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1130  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1131  if ( Normal.dot(v) <= 0 ) { return sd; }
1132  }
1133  }
1134  }
1135  }
1136  else
1137  {
1138  if (b>0) { sd = -b - std::sqrt(d); }
1139  else { sd = c/(-b+std::sqrt(d)); }
1140  zi = p.z() + sd*v.z() ;
1141  ri = rMinAv + zi*tanRMin ;
1142 
1143  if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1144  {
1145  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1146  { // seen on 64 bits systems. Split and recompute
1147  G4double fTerm = sd-std::fmod(sd,dRmax);
1148  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1149  }
1150  if ( !fPhiFullCone )
1151  {
1152  xi = p.x() + sd*v.x() ;
1153  yi = p.y() + sd*v.y() ;
1154  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1155 
1156  if (cosPsi >= cosHDPhiIT)
1157  {
1158  if ( sd > halfRadTolerance ) { snxt=sd; }
1159  else
1160  {
1161  // Calculate a normal vector in order to check Direction
1162 
1163  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1164  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1165  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1166  }
1167  }
1168  }
1169  else
1170  {
1171  if ( sd > halfRadTolerance ) { return sd; }
1172  else
1173  {
1174  // Calculate a normal vector in order to check Direction
1175 
1176  xi = p.x() + sd*v.x() ;
1177  yi = p.y() + sd*v.y() ;
1178  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1179  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1180  if ( Normal.dot(v) <= 0 ) { return sd; }
1181  }
1182  }
1183  }
1184  }
1185  }
1186  }
1187  else
1188  {
1189  // Within kRadTol*0.5 of inner cone (real OR imaginary)
1190  // ----> Check not travelling through (=>0 to in)
1191  // ----> if not:
1192  // -2nd root with validity check
1193 
1194  if ( std::fabs(p.z()) <= tolODz )
1195  {
1196  if ( nt2 > 0 )
1197  {
1198  // Inside inner real cone, heading outwards, inside z range
1199 
1200  if ( !fPhiFullCone )
1201  {
1202  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1203 
1204  if (cosPsi >= cosHDPhiIT) { return 0.0; }
1205  }
1206  else { return 0.0; }
1207  }
1208  else
1209  {
1210  // Within z extent, but not travelling through
1211  // -> 2nd root or kInfinity if 1st root on imaginary cone
1212 
1213  b = nt2/nt1 ;
1214  c = nt3/nt1 ;
1215  d = b*b - c ;
1216 
1217  if ( d >= 0 ) // > 0
1218  {
1219  if (b>0) { sd = -b - std::sqrt(d); }
1220  else { sd = c/(-b+std::sqrt(d)); }
1221  zi = p.z() + sd*v.z() ;
1222  ri = rMinAv + zi*tanRMin ;
1223 
1224  if ( ri > 0 ) // 2nd root
1225  {
1226  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1227  else { sd = -b + std::sqrt(d); }
1228 
1229  zi = p.z() + sd*v.z() ;
1230 
1231  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1232  {
1233  if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1234  { // seen on 64 bits systems. Split and recompute
1235  G4double fTerm = sd-std::fmod(sd,dRmax);
1236  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1237  }
1238  if ( !fPhiFullCone )
1239  {
1240  xi = p.x() + sd*v.x() ;
1241  yi = p.y() + sd*v.y() ;
1242  ri = rMinAv + zi*tanRMin ;
1243  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1244 
1245  if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1246  }
1247  else { return sd; }
1248  }
1249  }
1250  else { return kInfinity; }
1251  }
1252  }
1253  }
1254  else // 2nd root
1255  {
1256  b = nt2/nt1 ;
1257  c = nt3/nt1 ;
1258  d = b*b - c ;
1259 
1260  if ( d > 0 )
1261  {
1262  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1263  else { sd = -b + std::sqrt(d) ; }
1264  zi = p.z() + sd*v.z() ;
1265 
1266  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1267  {
1268  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1269  { // seen on 64 bits systems. Split and recompute
1270  G4double fTerm = sd-std::fmod(sd,dRmax);
1271  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1272  }
1273  if ( !fPhiFullCone )
1274  {
1275  xi = p.x() + sd*v.x();
1276  yi = p.y() + sd*v.y();
1277  ri = rMinAv + zi*tanRMin ;
1278  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1279 
1280  if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1281  }
1282  else { return sd; }
1283  }
1284  }
1285  }
1286  }
1287  }
1288  }
1289 
1290  // Phi segment intersection
1291  //
1292  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1293  //
1294  // o NOTE: Large duplication of code between sphi & ephi checks
1295  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1296  // intersection check <=0 -> >=0
1297  // -> Should use some form of loop Construct
1298 
1299  if ( !fPhiFullCone )
1300  {
1301  // First phi surface (starting phi)
1302 
1303  Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1304 
1305  if ( Comp < 0 ) // Component in outwards normal dirn
1306  {
1307  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1308 
1309  if (Dist < halfCarTolerance)
1310  {
1311  sd = Dist/Comp ;
1312 
1313  if ( sd < snxt )
1314  {
1315  if ( sd < 0 ) { sd = 0.0; }
1316 
1317  zi = p.z() + sd*v.z() ;
1318 
1319  if ( std::fabs(zi) <= tolODz )
1320  {
1321  xi = p.x() + sd*v.x() ;
1322  yi = p.y() + sd*v.y() ;
1323  rhoi2 = xi*xi + yi*yi ;
1324  tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1325  tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1326 
1327  if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1328  {
1329  // z and r intersections good - check intersecting with
1330  // correct half-plane
1331 
1332  if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1333  }
1334  }
1335  }
1336  }
1337  }
1338 
1339  // Second phi surface (Ending phi)
1340 
1341  Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1342 
1343  if ( Comp < 0 ) // Component in outwards normal dirn
1344  {
1345  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1346  if (Dist < halfCarTolerance)
1347  {
1348  sd = Dist/Comp ;
1349 
1350  if ( sd < snxt )
1351  {
1352  if ( sd < 0 ) { sd = 0.0; }
1353 
1354  zi = p.z() + sd*v.z() ;
1355 
1356  if (std::fabs(zi) <= tolODz)
1357  {
1358  xi = p.x() + sd*v.x() ;
1359  yi = p.y() + sd*v.y() ;
1360  rhoi2 = xi*xi + yi*yi ;
1361  tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1362  tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1363 
1364  if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1365  {
1366  // z and r intersections good - check intersecting with
1367  // correct half-plane
1368 
1369  if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1370  }
1371  }
1372  }
1373  }
1374  }
1375  }
1376  if (snxt < halfCarTolerance) { snxt = 0.; }
1377 
1378  return snxt ;
1379 }
1380 
1381 //////////////////////////////////////////////////////////////////////////////
1382 //
1383 // Calculate distance (<= actual) to closest surface of shape from outside
1384 // - Calculate distance to z, radial planes
1385 // - Only to phi planes if outside phi extent
1386 // - Return 0 if point inside
1387 
1388 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
1389 {
1390  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1391  G4double tanRMin, secRMin, pRMin ;
1392  G4double tanRMax, secRMax, pRMax ;
1393 
1394  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1395  safeZ = std::fabs(p.z()) - fDz ;
1396 
1397  if ( fRmin1 || fRmin2 )
1398  {
1399  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1400  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1401  pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1402  safeR1 = (pRMin - rho)/secRMin ;
1403 
1404  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1405  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1406  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1407  safeR2 = (rho - pRMax)/secRMax ;
1408 
1409  if ( safeR1 > safeR2) { safe = safeR1; }
1410  else { safe = safeR2; }
1411  }
1412  else
1413  {
1414  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1415  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1416  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1417  safe = (rho - pRMax)/secRMax ;
1418  }
1419  if ( safeZ > safe ) { safe = safeZ; }
1420 
1421  if ( !fPhiFullCone && rho )
1422  {
1423  // Psi=angle from central phi to point
1424 
1425  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1426 
1427  if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1428  {
1429  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1430  {
1431  safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1432  }
1433  else
1434  {
1435  safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1436  }
1437  if ( safePhi > safe ) { safe = safePhi; }
1438  }
1439  }
1440  if ( safe < 0.0 ) { safe = 0.0; }
1441 
1442  return safe ;
1443 }
1444 
1445 ///////////////////////////////////////////////////////////////
1446 //
1447 // Calculate distance to surface of shape from 'inside', allowing for tolerance
1448 // - Only Calc rmax intersection if no valid rmin intersection
1449 
1450 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
1451  const G4ThreeVector& v,
1452  const G4bool calcNorm,
1453  G4bool *validNorm,
1454  G4ThreeVector *n) const
1455 {
1456  ESide side = kNull, sider = kNull, sidephi = kNull;
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.y() != 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 
2069 //////////////////////////////////////////////////////////////////
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 
2147 ////////////////////////////////////////////////////////////////////////////
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 
2241 //////////////////////////////////////////////////////////////////////////
2242 //
2243 // GetEntityType
2244 
2245 G4GeometryType G4Cons::GetEntityType() const
2246 {
2247  return G4String("G4Cons");
2248 }
2249 
2250 //////////////////////////////////////////////////////////////////////////
2251 //
2252 // Make a clone of the object
2253 //
2254 G4VSolid* G4Cons::Clone() const
2255 {
2256  return new G4Cons(*this);
2257 }
2258 
2259 //////////////////////////////////////////////////////////////////////////
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 
2286 /////////////////////////////////////////////////////////////////////////
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 
2373 //////////////////////////////////////////////////////////////////////////
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 #endif
G4Cons::G4Cons
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Cons.cc:79
G4VSolid::GetName
G4String GetName() const
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4Polyhedron.hh
kMZ
Definition: G4Cons.cc:68
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VGraphicsScene.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
kRMin
Definition: G4Cons.cc:68
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4Cons::operator=
G4Cons & operator=(const G4Cons &rhs)
Definition: G4Cons.cc:170
G4VoxelLimits.hh
G4Cons::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Cons.cc:718
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:203
p
const char * p
Definition: xmltok.h:285
G4Cons::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Cons.cc:1450
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4Cons::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Cons.cc:2263
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
kMinMeshSections
const G4int kMinMeshSections
Definition: meshdefs.hh:45
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4Cons::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Cons.cc:2245
G4Cons::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Cons.cc:258
CLHEP::Hep3Vector::z
double z() const
kNull
Definition: G4Cons.cc:68
G4VSolid::DumpInfo
void DumpInfo() const
ENorm
ENorm
Definition: G4Cons.cc:72
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4cout
G4GLOB_DLL std::ostream G4cout
kSPhi
Definition: G4Cons.cc:68
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4bool
bool G4bool
Definition: G4Types.hh:79
G4Cons
Definition: G4Cons.hh:82
G4VPVParameterisation.hh
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
Randomize.hh
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4AffineTransform.hh
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4Cons::~G4Cons
~G4Cons()
Definition: G4Cons.cc:144
EInside
EInside
Definition: geomdefs.hh:58
EAxis
EAxis
Definition: geomdefs.hh:54
G4Cons::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Cons.cc:2382
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
plottest35.t1
tuple t1
Definition: plottest35.py:33
CLHEP::Hep3Vector::y
double y() const
FatalException
Definition: G4ExceptionSeverity.hh:60
kPZ
Definition: G4Cons.cc:68
G4Cons::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Cons.cc:2377
G4Cons.hh
G4Cons::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Cons.cc:270
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:305
kNEPhi
Definition: G4Cons.cc:72
ESide
ESide
Definition: G4Cons.cc:68
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4Cons::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Cons.cc:2290
kRMax
Definition: G4Cons.cc:68
G4Cons::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Cons.cc:488
G4double
double G4double
Definition: G4Types.hh:76
G4Cons::Clone
G4VSolid * Clone() const
Definition: G4Cons.cc:2254
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
G4GeometryTolerance.hh
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
kNRMax
Definition: G4Cons.cc:72
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4VSolid::ClipBetweenSections
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
kEPhi
Definition: G4Cons.cc:68
kMeshAngleDefault
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:74
kNRMin
Definition: G4Cons.cc:72
kMaxMeshSections
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
kNSPhi
Definition: G4Cons.cc:72
kNZ
Definition: G4Cons.cc:72
Generated on Wed Apr 30 2014 15:54:40 for Geant4.10 by   doxygen 1.8.7