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