source: trunk/source/geometry/solids/CSG/src/G4Cons.cc @ 916

Last change on this file since 916 was 850, checked in by garnier, 16 years ago

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