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

Last change on this file was 1337, checked in by garnier, 14 years ago

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