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