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

Last change on this file since 921 was 921, checked in by garnier, 15 years ago

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