source: trunk/source/geometry/solids/specific/src/G4GenericTrap.cc @ 1316

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 55.0 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: G4GenericTrap.cc,v 1.12 2010/06/11 09:42:28 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34// G4GenericTrap.cc
35//
36// Authors:
37//   Tatiana Nikitina, CERN; Ivana Hrivnacova, IPN Orsay
38//   Adapted from Root Arb8 implementation by Andrei Gheata, CERN
39// --------------------------------------------------------------------
40
41#include <iomanip>
42
43#include "G4GenericTrap.hh"
44#include "G4TessellatedSolid.hh"
45#include "G4TriangularFacet.hh"
46#include "G4QuadrangularFacet.hh"
47#include "Randomize.hh"
48
49#include "G4VGraphicsScene.hh"
50#include "G4Polyhedron.hh"
51#include "G4PolyhedronArbitrary.hh"
52#include "G4NURBS.hh"
53#include "G4NURBSbox.hh"
54#include "G4VisExtent.hh"
55
56const G4int    G4GenericTrap::fgkNofVertices = 8;
57const G4double G4GenericTrap::fgkTolerance = 1E-3;
58
59// --------------------------------------------------------------------
60
61G4GenericTrap::G4GenericTrap( const G4String& name, G4double hz,
62                              const std::vector<G4TwoVector>&  vertices )
63  : G4VSolid(name),
64    fpPolyhedron(0),
65    fDz(hz),
66    fVertices(),
67    fIsTwisted(false),
68    fTessellatedSolid(0),
69    fMinBBoxVector(G4ThreeVector(0,0,0)),
70    fMaxBBoxVector(G4ThreeVector(0,0,0)),
71    fVisSubdivisions(0),
72    fSurfaceArea(0.),
73    fCubicVolume(0.)
74   
75{
76  // General constructor
77
78  G4String errorDescription = "InvalidSetup in \" ";
79  errorDescription += name;
80  errorDescription += "\""; 
81 
82  // Check vertices size
83
84  if ( G4int(vertices.size()) != fgkNofVertices )
85  {
86    G4Exception("G4GenericTrap::G4GenericTrap()", errorDescription,
87                FatalException, "Number of vertices != 8");
88  }           
89   
90  // Check Ordering and Copy vertices
91  //
92  if(CheckOrder(vertices))
93  {
94    for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
95  }
96  else
97  { 
98    for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
99    for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
100  }
101
102  // Compute Twist
103  //
104  for( G4int i=0; i<4; i++) { fTwist[i]=0; }
105  fIsTwisted = ComputeIsTwisted();
106
107  // Compute Bounding Box
108  //
109  ComputeBBox();
110 
111  // If not twisted - create tessellated solid
112  // (an alternative implementation for testing)
113  //
114#ifdef G4TESS_TEST
115   if ( !fIsTwisted )  { fTessellatedSolid = CreateTessellatedSolid(); }
116#endif
117}
118
119// --------------------------------------------------------------------
120
121G4GenericTrap::G4GenericTrap( __void__& a )
122  : G4VSolid(a),
123    fpPolyhedron(0),
124    fDz(0.),
125    fVertices(),
126    fIsTwisted(false),
127    fTessellatedSolid(0),
128    fMinBBoxVector(G4ThreeVector(0,0,0)),
129    fMaxBBoxVector(G4ThreeVector(0,0,0)),
130    fVisSubdivisions(0),
131    fSurfaceArea(0.),
132    fCubicVolume(0.)
133   
134{
135  // Fake default constructor - sets only member data and allocates memory
136  //                            for usage restricted to object persistency.
137}
138
139// --------------------------------------------------------------------
140
141G4GenericTrap::~G4GenericTrap()
142{
143  // Destructor
144}
145
146// --------------------------------------------------------------------
147
148EInside
149G4GenericTrap::InsidePolygone(const G4ThreeVector& p,
150                              const std::vector<G4TwoVector>& poly) const 
151{
152  static const G4double halfCarTolerance=kCarTolerance*0.5;
153  EInside  in = kInside;
154  G4double cross, len2;
155  G4int count=0;
156
157  for (G4int i = 0; i < 4; i++)
158  {
159    G4int j = (i+1) % 4;
160
161    cross = (p.x()-poly[i].x())*(poly[j].y()-poly[i].y())-
162            (p.y()-poly[i].y())*(poly[j].x()-poly[i].x());
163
164    len2=(poly[i]-poly[j]).mag2();
165    if (len2 > kCarTolerance)
166    {
167      if(cross*cross<=len2*halfCarTolerance*halfCarTolerance)  // Surface check
168      {
169        G4double test;
170        G4int k,l;
171        if ( poly[i].y() > poly[j].y() )  { k=i; l=j; }
172        else                              { k=j; l=i; }
173
174        if ( poly[k].x() != poly[l].x() )
175        {
176          test = (p.x()-poly[l].x())/(poly[k].x()-poly[l].x())
177               * (poly[k].y()-poly[l].y())+poly[l].y();
178        }
179        else
180        {
181          test = p.y();
182        }
183
184        // Check if point is Inside Segment
185        //
186        if( (test>=(poly[l].y()-halfCarTolerance))
187         && (test<=(poly[k].y()+halfCarTolerance)) )
188        { 
189          return kSurface;
190        }
191        else
192        {
193          return kOutside;
194        }
195      }
196      else if (cross<0.)  { return kOutside; }
197    }
198    else
199    {
200      count++;
201    }
202  }
203
204  // All collapsed vertices, Tet like
205  //
206  if(count==4)
207  { 
208    if ( (fabs(p.x()-poly[0].x())+fabs(p.y()-poly[0].y())) > halfCarTolerance )
209    {
210      in=kOutside;
211    }
212  }
213  return in;
214}
215
216// --------------------------------------------------------------------
217
218EInside G4GenericTrap::Inside(const G4ThreeVector& p) const
219{
220  // Test if point is inside this shape
221
222#ifdef G4TESS_TEST
223   if ( fTessellatedSolid )
224   { 
225     return fTessellatedSolid->Inside(p);
226   }
227#endif 
228
229  static const G4double halfCarTolerance=kCarTolerance*0.5;
230  EInside innew=kOutside;
231  std::vector<G4TwoVector> xy;
232 
233  if (std::fabs(p.z()) <= fDz+halfCarTolerance)  // First check Z range
234  {
235    // Compute intersection between Z plane containing point and the shape
236    //
237    G4double cf = 0.5*(fDz-p.z())/fDz;
238    for (G4int i=0; i<4; i++)
239    {
240      xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
241    }
242
243    innew=InsidePolygone(p,xy);
244
245    if( (innew==kInside) || (innew==kSurface) )
246    { 
247      if(std::fabs(p.z()) > fDz-halfCarTolerance)  { innew=kSurface; }
248    }
249  }
250  return innew;   
251} 
252
253// --------------------------------------------------------------------
254
255G4ThreeVector G4GenericTrap::SurfaceNormal( const G4ThreeVector& p ) const
256{
257  // Calculate side nearest to p, and return normal
258  // If two sides are equidistant, sum of the Normal is returned
259
260#ifdef G4TESS_TEST
261  if ( fTessellatedSolid )
262  {
263    return fTessellatedSolid->SurfaceNormal(p);
264  } 
265#endif   
266
267  static const G4double halfCarTolerance=kCarTolerance*0.5;
268 
269  G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
270                p0, p1, p2, r1, r2, r3, r4;
271  G4int noSurfaces = 0; 
272  G4double distxy,distz;
273  G4bool zPlusSide=false;
274
275  distz = fDz-std::fabs(p.z());
276  if (distz < halfCarTolerance)
277  {
278    if(p.z()>0)
279    {
280      zPlusSide=true;
281      sumnorm=G4ThreeVector(0,0,1);
282    }
283    else
284    {
285      sumnorm=G4ThreeVector(0,0,-1);
286    }
287    noSurfaces ++;
288  } 
289
290  // Check lateral planes
291  //
292  std:: vector<G4TwoVector> vertices; 
293  G4double cf = 0.5*(fDz-p.z())/fDz;
294  for (G4int i=0; i<4; i++)
295  {
296    vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
297  }
298
299  // Compute distance for lateral planes
300  //
301  for (G4int s=0; s<4; s++)
302  {
303    p0=G4ThreeVector(vertices[s].x(),vertices[s].y(),p.z());
304    if(zPlusSide)
305    {
306      p1=G4ThreeVector(fVertices[s].x(),fVertices[s].y(),-fDz);
307    }
308    else
309    {
310      p1=G4ThreeVector(fVertices[s+4].x(),fVertices[s+4].y(),fDz); 
311    }
312    p2=G4ThreeVector(vertices[(s+1)%4].x(),vertices[(s+1)%4].y(),p.z());
313
314    // Collapsed vertices
315    //
316    if ( (p2-p0).mag2() < kCarTolerance )
317    {
318      if ( fabs(p.z()+fDz) > kCarTolerance )
319      {
320        p2=G4ThreeVector(fVertices[(s+1)%4].x(),fVertices[(s+1)%4].y(),-fDz);
321      }
322      else
323      {
324        p2=G4ThreeVector(fVertices[(s+1)%4+4].x(),fVertices[(s+1)%4+4].y(),fDz);
325      }
326    }
327    lnorm = (p1-p0).cross(p2-p0);
328    lnorm = lnorm.unit();
329
330    // Adjust Normal for Twisted Surface
331    //
332    if ( (fIsTwisted) && (GetTwistAngle(s)!=0) )
333    {
334      G4double normP=(p2-p0).mag();
335      if(normP)
336      {
337        G4double proj=(p-p0).dot(p2-p0)/normP;
338        if(proj<0)     { proj=0; }
339        if(proj>normP) { proj=normP; }
340        G4int j=(s+1)%4;
341        r1=G4ThreeVector(fVertices[s+4].x(),fVertices[s+4].y(),fDz);
342        r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
343        r3=G4ThreeVector(fVertices[s].x(),fVertices[s].y(),-fDz);
344        r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
345        r1=r1+proj*(r2-r1)/normP;
346        r3=r3+proj*(r4-r3)/normP;
347        r2=r1-r3;
348        r4=r2.cross(p2-p0); r4=r4.unit();
349        lnorm=r4;
350      }
351    }   // End if fIsTwisted
352
353    distxy=std::fabs((p0-p).dot(lnorm));
354    if ( distxy<halfCarTolerance )
355    {
356      noSurfaces ++;
357
358      // Negative sign for Normal is taken for Outside Normal
359      //
360      sumnorm=sumnorm+lnorm;
361    }
362
363    // For ApproxSurfaceNormal
364    //
365    if (distxy<distz)
366    {
367      distz=distxy;
368      apprnorm=lnorm;
369    }     
370  }  // End for loop
371
372  // Calculate final Normal, add Normal in the Corners and Touching Sides
373  //
374  if ( noSurfaces == 0 )
375  {
376    G4Exception("G4GenericTrap::SurfaceNormal(p)", "Notification",
377                JustWarning, "Point p is not on surface !?" );
378    sumnorm=apprnorm;
379    // Add Approximative Surface Normal Calculation?
380  }
381  else if ( noSurfaces == 1 )  { sumnorm = sumnorm; }
382  else                         { sumnorm = sumnorm.unit(); }
383
384  return sumnorm ; 
385}   
386
387// --------------------------------------------------------------------
388
389G4ThreeVector G4GenericTrap::NormalToPlane( const G4ThreeVector& p,
390                                            const G4int ipl ) const
391{
392  // Return normal to given lateral plane ipl
393
394#ifdef G4TESS_TEST
395  if ( fTessellatedSolid )
396  { 
397    return fTessellatedSolid->SurfaceNormal(p);
398  } 
399#endif   
400
401  static const G4double halfCarTolerance=kCarTolerance*0.5; 
402  G4ThreeVector lnorm, norm(0.,0.,0.), p0,p1,p2;
403 
404  G4double  distz = fDz-p.z();
405  G4int i=ipl;  // current plane index
406 
407  G4TwoVector u,v; 
408  G4ThreeVector r1,r2,r3,r4;
409  G4double cf = 0.5*(fDz-p.z())/fDz;
410  G4int j=(i+1)%4;
411
412  u=fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]);
413  v=fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
414
415  // Compute cross product
416  //
417  p0=G4ThreeVector(u.x(),u.y(),p.z());
418     
419  if (std::fabs(distz)<halfCarTolerance)
420  {
421    p1=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);distz=-1;}
422  else
423  {
424    p1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
425  } 
426  p2=G4ThreeVector(v.x(),v.y(),p.z());
427
428  // Collapsed vertices
429  //
430  if ( (p2-p0).mag2() < kCarTolerance )
431  {
432    if ( fabs(p.z()+fDz) > halfCarTolerance )
433    {
434      p2=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
435    }
436    else
437    {
438      p2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
439    }
440  }
441  lnorm=-(p1-p0).cross(p2-p0);
442  if (distz>-halfCarTolerance)  { lnorm=-lnorm.unit(); }
443  else                          { lnorm=lnorm.unit();  }
444 
445  // Adjust Normal for Twisted Surface
446  //
447  if( (fIsTwisted) && (GetTwistAngle(ipl)!=0) )
448  {
449    G4double normP=(p2-p0).mag();
450    if(normP)
451    {
452      G4double proj=(p-p0).dot(p2-p0)/normP;
453      if (proj<0)     { proj=0; }
454      if (proj>normP) { proj=normP; }
455
456      r1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
457      r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
458      r3=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);
459      r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
460      r1=r1+proj*(r2-r1)/normP;
461      r3=r3+proj*(r4-r3)/normP;
462      r2=r1-r3;
463      r4=r2.cross(p2-p0);r4=r4.unit();
464      lnorm=r4;
465    }
466  }  // End if fIsTwisted
467
468  return lnorm;
469}   
470
471// --------------------------------------------------------------------
472
473G4double G4GenericTrap::DistToPlane(const G4ThreeVector& p,
474                                    const G4ThreeVector& v,
475                                    const G4int ipl) const 
476{
477  // Computes distance to plane ipl :
478  // ipl=0 : points 0,4,1,5
479  // ipl=1 : points 1,5,2,6
480  // ipl=2 : points 2,6,3,7
481  // ipl=3 : points 3,7,0,4
482
483  static const G4double halfCarTolerance=0.5*kCarTolerance;
484  G4double xa,xb,xc,xd,ya,yb,yc,yd;
485 
486  G4int j = (ipl+1)%4;
487 
488  xa=fVertices[ipl].x();
489  ya=fVertices[ipl].y();
490  xb=fVertices[ipl+4].x();
491  yb=fVertices[ipl+4].y();
492  xc=fVertices[j].x();
493  yc=fVertices[j].y();
494  xd=fVertices[4+j].x();
495  yd=fVertices[4+j].y();
496 
497  G4double dz2 =0.5/fDz;
498  G4double tx1 =dz2*(xb-xa);
499  G4double ty1 =dz2*(yb-ya);
500  G4double tx2 =dz2*(xd-xc);
501  G4double ty2 =dz2*(yd-yc);
502  G4double dzp =fDz+p.z();
503  G4double xs1 =xa+tx1*dzp;
504  G4double ys1 =ya+ty1*dzp;
505  G4double xs2 =xc+tx2*dzp;
506  G4double ys2 =yc+ty2*dzp;
507  G4double dxs =xs2-xs1;
508  G4double dys =ys2-ys1;
509  G4double dtx =tx2-tx1;
510  G4double dty =ty2-ty1;
511
512  G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
513  G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
514             + tx1*ys2-tx2*ys1)*v.z();
515  G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
516  G4double s=kInfinity;
517  G4double x1,x2,y1,y2,xp,yp,zi;
518
519  if (std::fabs(a)<kCarTolerance)
520  {           
521    if (std::fabs(b)<kCarTolerance)  { return kInfinity; }
522    s=-c/b;
523
524    // Check if Point is on the Surface
525
526    if (s>-halfCarTolerance)
527    {
528      if (s<halfCarTolerance)
529      {
530        if (NormalToPlane(p,ipl).dot(v)<=0) { return 0.; }
531        else                                { return kInfinity; }
532      }
533
534      // Check the Intersection
535      //
536      zi=p.z()+s*v.z();
537      if (std::fabs(zi)<fDz)
538      {
539        x1=xs1+tx1*v.z()*s;
540        x2=xs2+tx2*v.z()*s;
541        xp=p.x()+s*v.x();
542        y1=ys1+ty1*v.z()*s;
543        y2=ys2+ty2*v.z()*s;
544        yp=p.y()+s*v.y();
545        zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
546        if (zi<=halfCarTolerance)  { return s; }
547      }
548    }
549    return kInfinity;
550  }     
551  G4double d=b*b-4*a*c;
552  if (d>=0)
553  {
554    if (a>0) { s=0.5*(-b-std::sqrt(d))/a; }
555    else     { s=0.5*(-b+std::sqrt(d))/a; }
556
557    // Check if Point is on the Surface
558    //
559    if (s>-halfCarTolerance)
560    {
561      if(s<halfCarTolerance)
562      {
563        if (NormalToPlane(p,ipl).dot(v)<=0)  { return 0.;}
564        else  // Check second root; return kInfinity
565        {
566          if (a>0) { s=0.5*(-b+std::sqrt(d))/a; }
567          else     { s=0.5*(-b-std::sqrt(d))/a; }
568          if (s<=halfCarTolerance) { return kInfinity; }
569        }
570      }
571      // Check the Intersection
572      //
573      zi=p.z()+s*v.z();
574      if (std::fabs(zi)<fDz)
575      {
576        x1=xs1+tx1*v.z()*s;
577        x2=xs2+tx2*v.z()*s;
578        xp=p.x()+s*v.x();
579        y1=ys1+ty1*v.z()*s;
580        y2=ys2+ty2*v.z()*s;
581        yp=p.y()+s*v.y();
582        zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
583        if (zi<=halfCarTolerance)  { return s; }
584      }
585    }
586    if (a>0)  { s=0.5*(-b+std::sqrt(d))/a; }
587    else      { s=0.5*(-b-std::sqrt(d))/a; }
588
589    // Check if Point is on the Surface
590    //
591    if (s>-halfCarTolerance)
592    {
593      if(s<halfCarTolerance)
594      {
595        if (NormalToPlane(p,ipl).dot(v)<=0)  { return 0.; }
596        else   // Check second root; return kInfinity.
597        {
598          if (a>0) { s=0.5*(-b-std::sqrt(d))/a; }
599          else     { s=0.5*(-b+std::sqrt(d))/a; }
600          if (s<=halfCarTolerance)  { return kInfinity; }
601        }
602      }
603      // Check the Intersection
604      //
605      zi=p.z()+s*v.z();
606      if (std::fabs(zi)<fDz)
607      {
608        x1=xs1+tx1*v.z()*s;
609        x2=xs2+tx2*v.z()*s;
610        xp=p.x()+s*v.x();
611        y1=ys1+ty1*v.z()*s;
612        y2=ys2+ty2*v.z()*s;
613        yp=p.y()+s*v.y();
614        zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
615        if (zi<=halfCarTolerance)  { return s; }
616      }
617    }
618  }
619  return kInfinity;
620}     
621
622// --------------------------------------------------------------------
623
624G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p,
625                                     const G4ThreeVector& v) const
626{
627#ifdef G4TESS_TEST
628  if ( fTessellatedSolid )
629  { 
630    return fTessellatedSolid->DistanceToIn(p, v);
631  } 
632#endif
633   
634  static const G4double halfCarTolerance=kCarTolerance*0.5;
635
636  G4double dist[5];
637  G4ThreeVector n;
638
639  // Check lateral faces
640  //
641  G4int i;
642  for (i=0; i<4; i++)
643  {
644    dist[i]=DistToPlane(p, v, i); 
645  }
646
647  // Check Z planes
648  //
649  dist[4]=kInfinity;
650  if (std::fabs(p.z())>fDz-halfCarTolerance)
651  {
652    if (v.z())
653    {
654      G4ThreeVector pt;
655      if (p.z()>0)
656      {
657        dist[4] = (fDz-p.z())/v.z();
658      }
659      else
660      {   
661        dist[4] = (-fDz-p.z())/v.z();
662      }   
663      if (dist[4]<-halfCarTolerance)
664      {
665        dist[4]=kInfinity;
666      }
667      else
668      {
669        if(dist[4]<halfCarTolerance)
670        {
671          if(p.z()>0)  { n=G4ThreeVector(0,0,1); }
672          else         { n=G4ThreeVector(0,0,-1); }
673          if (n.dot(v)<0) { dist[4]=0.; }
674          else            { dist[4]=kInfinity; }
675        }
676        pt=p+dist[4]*v;
677        if (Inside(pt)==kOutside)  { dist[4]=kInfinity; }
678      }
679    }
680  }   
681  G4double distmin = dist[0];
682  for (i=1;i<5;i++)
683  {
684    if (dist[i] < distmin)  { distmin = dist[i]; }
685  }
686
687  if (distmin<halfCarTolerance)  { distmin=0.; }
688
689  return distmin;
690}   
691 
692// --------------------------------------------------------------------
693
694G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p) const
695{
696  // Computes the closest distance from given point to this shape
697
698#ifdef G4TESS_TEST
699  if ( fTessellatedSolid )
700  {
701    return fTessellatedSolid->DistanceToIn(p);
702  } 
703#endif
704 
705  G4double safz = std::fabs(p.z())-fDz;
706  if(safz<0) { safz=0; }
707
708  G4int iseg;
709  G4double safe  = safz;
710  G4double safxy = safz;
711 
712  for (iseg=0; iseg<4; iseg++)
713  { 
714    safxy = SafetyToFace(p,iseg);
715    if (safxy>safe)  { safe=safxy; }
716  }
717
718  return safe;
719}
720
721// --------------------------------------------------------------------
722
723G4double
724G4GenericTrap::SafetyToFace(const G4ThreeVector& p, const G4int iseg) const
725{
726  // Estimate distance to lateral plane defined by segment iseg in range [0,3]
727  // Might be negative: plane seen only from inside
728
729  G4ThreeVector p1,norm;
730  G4double safe;
731
732  p1=G4ThreeVector(fVertices[iseg].x(),fVertices[iseg].y(),-fDz);
733 
734  norm=NormalToPlane(p,iseg);
735  safe = (p-p1).dot(norm); // Can be negative
736 
737  return safe;
738}
739
740// --------------------------------------------------------------------
741
742G4double
743G4GenericTrap::DistToTriangle(const G4ThreeVector& p,
744                              const G4ThreeVector& v, const G4int ipl) const
745{
746  static const G4double halfCarTolerance=kCarTolerance*0.5;
747
748  G4double xa=fVertices[ipl].x();
749  G4double ya=fVertices[ipl].y();
750  G4double xb=fVertices[ipl+4].x();
751  G4double yb=fVertices[ipl+4].y();
752  G4int j=(ipl+1)%4;
753  G4double xc=fVertices[j].x();
754  G4double yc=fVertices[j].y();
755  G4double zab=2*fDz;
756  G4double zac=0;
757 
758  if ( (fabs(xa-xc)+fabs(ya-yc)) < halfCarTolerance )
759  {
760    xc=fVertices[j+4].x();
761    yc=fVertices[j+4].y();
762    zac=2*fDz;
763    zab=2*fDz;
764
765    //Line case
766    //
767    if ( (fabs(xb-xc)+fabs(yb-yc)) < halfCarTolerance )
768    {
769      return kInfinity;
770    }
771  }
772  G4double a=(yb-ya)*zac-(yc-ya)*zab;
773  G4double b=(xc-xa)*zab-(xb-xa)*zac;
774  G4double c=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
775  G4double d=-xa*a-ya*b+fDz*c;
776  G4double t=a*v.x()+b*v.y()+c*v.z();
777
778  if (t!=0)
779  {
780    t=-(a*p.x()+b*p.y()+c*p.z()+d)/t;
781  }
782  if ( (t<halfCarTolerance) && (t>-halfCarTolerance) )
783  {
784    if (NormalToPlane(p,ipl).dot(v)<0)
785    {
786      t=kInfinity;
787    }
788    else
789    {
790      t=0;
791    }
792  }
793  if (Inside(p+v*t) != kSurface)  { t=kInfinity; }
794 
795  return t;
796} 
797
798// --------------------------------------------------------------------
799
800G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p,
801                                      const G4ThreeVector& v,
802                                      const G4bool calcNorm,
803                                            G4bool* validNorm,
804                                            G4ThreeVector* n) const
805{
806#ifdef G4TESS_TEST
807  if ( fTessellatedSolid )
808  { 
809    return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
810  }
811#endif
812
813  static const G4double halfCarTolerance=kCarTolerance*0.5;
814
815  G4double distmin;
816  G4bool lateral_cross = false;
817  ESide side = kUndefined;
818 
819  if (calcNorm)  { *validNorm=true; } // All normals are valid
820
821  if (v.z() < 0)
822  {
823    distmin=(-fDz-p.z())/v.z();
824    if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
825  }
826  else
827  {
828    if (v.z() > 0)
829    { 
830      distmin = (fDz-p.z())/v.z(); 
831      if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); } 
832    }
833    else            { distmin = kInfinity; }
834  }     
835
836  G4double dz2 =0.5/fDz;
837  G4double xa,xb,xc,xd;
838  G4double ya,yb,yc,yd;
839
840  for (G4int ipl=0; ipl<4; ipl++)
841  {
842    G4int j = (ipl+1)%4;
843    xa=fVertices[ipl].x();
844    ya=fVertices[ipl].y();
845    xb=fVertices[ipl+4].x();
846    yb=fVertices[ipl+4].y();
847    xc=fVertices[j].x();
848    yc=fVertices[j].y();
849    xd=fVertices[4+j].x();
850    yd=fVertices[4+j].y();
851
852    if ( ((fabs(xb-xd)+fabs(yb-yd))<halfCarTolerance)
853      || ((fabs(xa-xc)+fabs(ya-yc))<halfCarTolerance) )
854    {
855      G4double s=DistToTriangle(p,v,ipl) ;
856      if ( (s>=0) && (s<distmin) )
857      {
858        distmin=s;
859        lateral_cross=true;
860        side=ESide(ipl+1);
861      }
862      continue;
863    }
864    G4double tx1 =dz2*(xb-xa);
865    G4double ty1 =dz2*(yb-ya);
866    G4double tx2 =dz2*(xd-xc);
867    G4double ty2 =dz2*(yd-yc);
868    G4double dzp =fDz+p.z();
869    G4double xs1 =xa+tx1*dzp;
870    G4double ys1 =ya+ty1*dzp;
871    G4double xs2 =xc+tx2*dzp;
872    G4double ys2 =yc+ty2*dzp;
873    G4double dxs =xs2-xs1;
874    G4double dys =ys2-ys1;
875    G4double dtx =tx2-tx1;
876    G4double dty =ty2-ty1;
877    G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
878    G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
879               + tx1*ys2-tx2*ys1)*v.z();
880    G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
881    G4double s=kInfinity;
882
883    if (std::fabs(a) < kCarTolerance)
884    {           
885      if (std::fabs(b) < kCarTolerance) { continue; }
886      s=-c/b;
887         
888      // Check for Point on the Surface
889      //
890      if ((s > -halfCarTolerance) && (s < distmin))
891      {
892        if (s < halfCarTolerance)
893        {
894          if (NormalToPlane(p,ipl).dot(v)<0.)  { continue; }
895        }
896        distmin =s;
897        lateral_cross=true;
898        side=ESide(ipl+1);
899      }   
900      continue;
901    }
902    G4double d=b*b-4*a*c;
903    if (d >= 0.)
904    {
905      if (a > 0)  { s=0.5*(-b-std::sqrt(d))/a; }
906      else        { s=0.5*(-b+std::sqrt(d))/a; }
907
908      // Check for Point on the Surface
909      //
910      if (s > -halfCarTolerance )
911      {
912        if (s < distmin)
913        {
914          if(s < halfCarTolerance)
915          {
916            if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
917            {
918              if (a > 0)  { s=0.5*(-b+std::sqrt(d))/a; }
919              else        { s=0.5*(-b-std::sqrt(d))/a; }
920              if (( s > halfCarTolerance) && (s < distmin))
921              {
922                distmin=s;
923                lateral_cross = true;
924                side=ESide(ipl+1);
925              }
926              continue;
927            }
928          }
929          distmin = s;
930          lateral_cross = true;
931          side=ESide(ipl+1);
932        }
933      }
934      else
935      {
936        if (a > 0)  { s=0.5*(-b+std::sqrt(d))/a; }
937        else        { s=0.5*(-b-std::sqrt(d))/a; }
938
939        // Check for Point on the Surface
940        //
941        if ((s > -halfCarTolerance) && (s < distmin))
942        {
943          if (s < halfCarTolerance)
944          {
945            if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
946            {
947              if (a > 0)  { s=0.5*(-b-std::sqrt(d))/a; }
948              else        { s=0.5*(-b+std::sqrt(d))/a; }
949              if ( ( s > halfCarTolerance) && (s < distmin) )
950              {
951                distmin=s;
952                lateral_cross = true;
953                side=ESide(ipl+1);
954              }
955              continue;
956            } 
957          }
958          distmin =s;
959          lateral_cross = true;
960          side=ESide(ipl+1);
961        }   
962      }
963    }
964  }
965  if (!lateral_cross)  // Make sure that track crosses the top or bottom
966  {
967    if (distmin >= kInfinity)  { distmin=kCarTolerance; }
968    G4ThreeVector pt=p+distmin*v;
969   
970    // Check if propagated point is in the polygon
971    //
972    G4int i=0;
973    if (v.z()>0.) { i=4; }
974    std::vector<G4TwoVector> xy;
975    for ( G4int j=0; j<4; j++)  { xy.push_back(fVertices[i+j]); }
976
977    // Check Inside
978    //
979    if (InsidePolygone(pt,xy)==kOutside)
980    { 
981      if(calcNorm)
982      {
983        if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
984        else         { side=kMZ; *n = G4ThreeVector(0,0,-1);}
985      } 
986      return 0.;
987    }
988    else
989    {
990      if(v.z()>0) {side=kPZ;}
991      else        {side=kMZ;}
992    }
993  }
994
995  if (calcNorm)
996  {
997    G4ThreeVector pt=p+v*distmin;     
998    switch (side)
999    {
1000      case kXY0:
1001        *n=NormalToPlane(pt,0);
1002        break;
1003      case kXY1:
1004        *n=NormalToPlane(pt,1);
1005        break;
1006      case kXY2:
1007        *n=NormalToPlane(pt,2);
1008        break;
1009      case kXY3:
1010        *n=NormalToPlane(pt,3);
1011        break;
1012      case kPZ:
1013        *n=G4ThreeVector(0,0,1);
1014        break;
1015      case kMZ:
1016        *n=G4ThreeVector(0,0,-1);
1017        break;
1018      default:
1019        G4cout.precision(16);
1020        G4cout << G4endl;
1021        DumpInfo();
1022        G4cout << "Position:"  << G4endl << G4endl;
1023        G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
1024        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
1025        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
1026        G4cout << "Direction:" << G4endl << G4endl;
1027        G4cout << "v.x() = "   << v.x() << G4endl;
1028        G4cout << "v.y() = "   << v.y() << G4endl;
1029        G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
1030        G4cout << "Proposed distance :" << G4endl << G4endl;
1031        G4cout << "distmin = "    << distmin/mm << " mm" << G4endl << G4endl;
1032        G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
1033                    "Notification", JustWarning,
1034                    "Undefined side for valid surface normal to solid.");
1035        break;
1036     }
1037  }
1038 
1039  if (distmin<halfCarTolerance)  { distmin=0.; }
1040 
1041  return distmin;
1042}   
1043
1044// --------------------------------------------------------------------
1045
1046G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p) const
1047{
1048
1049#ifdef G4TESS_TEST
1050  if ( fTessellatedSolid )
1051  { 
1052    return fTessellatedSolid->DistanceToOut(p);
1053  } 
1054#endif
1055
1056  G4double safz = fDz-std::fabs(p.z());
1057  if (safz<0) { safz = 0; }
1058
1059  G4double safe  = safz;
1060  G4double safxy = safz;
1061 
1062  for (G4int iseg=0; iseg<4; iseg++)
1063  { 
1064    safxy = std::fabs(SafetyToFace(p,iseg));
1065    if (safxy < safe)  { safe = safxy; }
1066  }
1067
1068  return safe;
1069}   
1070
1071// --------------------------------------------------------------------
1072
1073G4bool G4GenericTrap::CalculateExtent(const EAxis pAxis,
1074                                      const G4VoxelLimits& pVoxelLimit,
1075                                      const G4AffineTransform& pTransform,
1076                                      G4double& pMin, G4double& pMax) const
1077{
1078#ifdef G4TESS_TEST
1079  if ( fTessellatedSolid )
1080  {
1081    return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
1082                                              pTransform, pMin, pMax);
1083  }     
1084#endif
1085
1086  // Computes bounding vectors for a shape
1087  //
1088  G4double Dx,Dy;
1089  G4ThreeVector minVec = GetMinimumBBox();
1090  G4ThreeVector maxVec = GetMaximumBBox();
1091  Dx = 0.5*(maxVec.x()- minVec.y());
1092  Dy = 0.5*(maxVec.y()- minVec.y());
1093
1094  if (!pTransform.IsRotated())
1095  {
1096    // Special case handling for unrotated shapes
1097    // Compute x/y/z mins and maxs respecting limits, with early returns
1098    // if outside limits. Then switch() on pAxis
1099    //
1100    G4double xoffset,xMin,xMax;
1101    G4double yoffset,yMin,yMax;
1102    G4double zoffset,zMin,zMax;
1103
1104    xoffset=pTransform.NetTranslation().x();
1105    xMin=xoffset-Dx;
1106    xMax=xoffset+Dx;
1107    if (pVoxelLimit.IsXLimited())
1108    {
1109      if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
1110        || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
1111      {
1112        return false;
1113      }
1114      else
1115      {
1116        if (xMin<pVoxelLimit.GetMinXExtent())
1117        {
1118          xMin=pVoxelLimit.GetMinXExtent();
1119        }
1120        if (xMax>pVoxelLimit.GetMaxXExtent())
1121        {
1122          xMax=pVoxelLimit.GetMaxXExtent();
1123        }
1124      }
1125    }
1126
1127    yoffset=pTransform.NetTranslation().y();
1128    yMin=yoffset-Dy;
1129    yMax=yoffset+Dy;
1130    if (pVoxelLimit.IsYLimited())
1131    {
1132      if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
1133        || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
1134      {
1135        return false;
1136      }
1137      else
1138      {
1139        if (yMin<pVoxelLimit.GetMinYExtent())
1140        {
1141          yMin=pVoxelLimit.GetMinYExtent();
1142        }
1143        if (yMax>pVoxelLimit.GetMaxYExtent())
1144        {
1145          yMax=pVoxelLimit.GetMaxYExtent();
1146        }
1147      }
1148    }
1149
1150    zoffset=pTransform.NetTranslation().z();
1151    zMin=zoffset-fDz;
1152    zMax=zoffset+fDz;
1153    if (pVoxelLimit.IsZLimited())
1154    {
1155      if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
1156        || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
1157      {
1158        return false;
1159      }
1160      else
1161      {
1162        if (zMin<pVoxelLimit.GetMinZExtent())
1163        {
1164          zMin=pVoxelLimit.GetMinZExtent();
1165        }
1166        if (zMax>pVoxelLimit.GetMaxZExtent())
1167        {
1168          zMax=pVoxelLimit.GetMaxZExtent();
1169        }
1170      }
1171    }
1172
1173    switch (pAxis)
1174    {
1175      case kXAxis:
1176        pMin = xMin;
1177        pMax = xMax;
1178        break;
1179      case kYAxis:
1180        pMin = yMin;
1181        pMax = yMax;
1182        break;
1183      case kZAxis:
1184        pMin = zMin;
1185        pMax = zMax;
1186        break;
1187      default:
1188        break;
1189    }
1190    pMin-=kCarTolerance;
1191    pMax+=kCarTolerance;
1192
1193    return true;
1194  }
1195  else
1196  {
1197    // General rotated case - create and clip mesh to boundaries
1198
1199    G4bool existsAfterClip=false;
1200    G4ThreeVectorList *vertices;
1201
1202    pMin=+kInfinity;
1203    pMax=-kInfinity;
1204
1205    // Calculate rotated vertex coordinates
1206    //
1207    vertices=CreateRotatedVertices(pTransform);
1208    ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
1209    ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
1210    ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
1211
1212    if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
1213    {
1214      existsAfterClip=true;
1215   
1216      // Add 2*tolerance to avoid precision troubles
1217      //
1218      pMin-=kCarTolerance;
1219      pMax+=kCarTolerance;
1220    }
1221    else
1222    {
1223      // Check for case where completely enveloping clipping volume.
1224      // If point inside then we are confident that the solid completely
1225      // envelopes the clipping volume. Hence set min/max extents according
1226      // to clipping volume extents along the specified axis.
1227      //
1228      G4ThreeVector clipCentre(
1229              (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
1230              (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
1231              (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
1232
1233      if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
1234      {
1235        existsAfterClip=true;
1236        pMin=pVoxelLimit.GetMinExtent(pAxis);
1237        pMax=pVoxelLimit.GetMaxExtent(pAxis);
1238      }
1239    }
1240    delete vertices;
1241    return existsAfterClip;
1242  }
1243}
1244
1245// --------------------------------------------------------------------
1246
1247G4ThreeVectorList*
1248G4GenericTrap::CreateRotatedVertices(const G4AffineTransform& pTransform) const
1249{
1250  // Create a List containing the transformed vertices
1251  // Ordering [0-3] -fDz cross section
1252  //          [4-7] +fDz cross section such that [0] is below [4],
1253  //                                             [1] below [5] etc.
1254  // Note: caller has deletion responsibility
1255
1256  G4ThreeVector Min = GetMinimumBBox();
1257  G4ThreeVector Max = GetMaximumBBox();
1258
1259  G4ThreeVectorList *vertices;
1260  vertices=new G4ThreeVectorList();
1261  vertices->reserve(8);
1262   
1263  if (vertices)
1264  {
1265    G4ThreeVector vertex0(Min.x(),Min.y(),Min.z());
1266    G4ThreeVector vertex1(Max.x(),Min.y(),Min.z());
1267    G4ThreeVector vertex2(Max.x(),Max.y(),Min.z());
1268    G4ThreeVector vertex3(Min.x(),Max.y(),Min.z());
1269    G4ThreeVector vertex4(Min.x(),Min.y(),Max.z());
1270    G4ThreeVector vertex5(Max.x(),Min.y(),Max.z());
1271    G4ThreeVector vertex6(Max.x(),Max.y(),Max.z());
1272    G4ThreeVector vertex7(Min.x(),Max.y(),Max.z());
1273
1274    vertices->push_back(pTransform.TransformPoint(vertex0));
1275    vertices->push_back(pTransform.TransformPoint(vertex1));
1276    vertices->push_back(pTransform.TransformPoint(vertex2));
1277    vertices->push_back(pTransform.TransformPoint(vertex3));
1278    vertices->push_back(pTransform.TransformPoint(vertex4));
1279    vertices->push_back(pTransform.TransformPoint(vertex5));
1280    vertices->push_back(pTransform.TransformPoint(vertex6));
1281    vertices->push_back(pTransform.TransformPoint(vertex7));
1282  }
1283  else
1284  {
1285    G4Exception("G4GenericTrap::CreateRotatedVertices()", "FatalError",
1286                FatalException, "Out of memory - Cannot allocate vertices!");
1287  }
1288  return vertices;
1289}
1290 
1291// --------------------------------------------------------------------
1292
1293std::ostream& G4GenericTrap::StreamInfo(std::ostream& os) const
1294{
1295  os << "-----------------------------------------------------------\n"
1296     << "    *** Dump for solid - " << GetName() << " *** \n"
1297     << "    =================================================== \n"
1298     << " Solid geometry type: " << GetEntityType()  << G4endl
1299     << "   half length Z: " << fDz/mm << " mm \n"
1300     << "   list of vertices:\n";
1301     
1302  for ( G4int i=0; i<fgkNofVertices; ++i )
1303  {
1304    os << std::setw(5) << "#" << i
1305       << "   vx = " << fVertices[i].x()/mm << " mm" 
1306       << "   vy = " << fVertices[i].y()/mm << " mm" << G4endl;
1307  }
1308 
1309  return os;
1310} 
1311
1312// --------------------------------------------------------------------
1313
1314G4double G4GenericTrap::GetSurfaceArea()
1315{
1316  if(fSurfaceArea != 0.) {;}
1317  else
1318  {
1319    std::vector<G4ThreeVector> vertices;
1320    for (G4int i=0; i<4;i++)
1321    {
1322      vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1323    }
1324    for (G4int i=4; i<8;i++)
1325    {
1326      vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1327    }
1328
1329    // Surface Area of Planes(only estimation for twisted)
1330    //
1331    G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1332                                          vertices[2],vertices[3]);//-fDz plane
1333    G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1334                                          vertices[5],vertices[4]);// Lat plane
1335    G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1336                                          vertices[4],vertices[7]);// Lat plane
1337    G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1338                                          vertices[7],vertices[6]);// Lat plane
1339    G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1340                                          vertices[5],vertices[6]);// Lat plane
1341    G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1342                                          vertices[6],vertices[7]);// fDz plane
1343
1344    // Total Surface Area
1345    //
1346    if(!fIsTwisted)
1347    {
1348      fSurfaceArea = fSurface0+fSurface1+fSurface2
1349                   + fSurface3+fSurface4+fSurface5;
1350    }
1351    else
1352    {
1353      fSurfaceArea = G4VSolid::GetSurfaceArea();
1354    }
1355  }
1356  return fSurfaceArea;
1357}
1358
1359// --------------------------------------------------------------------
1360
1361G4double G4GenericTrap::GetFaceSurfaceArea(const G4ThreeVector& p0,
1362                                           const G4ThreeVector& p1, 
1363                                           const G4ThreeVector& p2,
1364                                           const G4ThreeVector& p3) const
1365{
1366  // Auxiliary method for Get Surface Area of Face
1367 
1368  G4double aOne, aTwo;
1369  G4ThreeVector t, u, v, w, Area, normal;
1370
1371  t = p2 - p1;
1372  u = p0 - p1;
1373  v = p2 - p3;
1374  w = p0 - p3;
1375 
1376  Area = w.cross(v);
1377  aOne = 0.5*Area.mag();
1378 
1379  Area = t.cross(u);
1380  aTwo = 0.5*Area.mag();
1381 
1382  return aOne + aTwo;
1383}
1384
1385// --------------------------------------------------------------------
1386
1387G4ThreeVector G4GenericTrap::GetPointOnSurface() const
1388{
1389
1390#ifdef G4TESS_TEST
1391  if ( fTessellatedSolid )
1392  { 
1393    return fTessellatedSolid->GetPointOnSurface();
1394  } 
1395#endif
1396
1397  G4ThreeVector point;
1398  G4TwoVector u,v,w;
1399  G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
1400  G4int ipl,j;
1401 
1402  std::vector<G4ThreeVector> vertices;
1403  for (G4int i=0; i<4;i++)
1404  {
1405    vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1406  }
1407  for (G4int i=4; i<8;i++)
1408  {
1409    vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1410  }
1411
1412  // Surface Area of Planes(only estimation for twisted)
1413  //
1414  G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1415                                       vertices[2],vertices[3]);//-fDz plane
1416  G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1417                                       vertices[5],vertices[4]);// Lat plane
1418  G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1419                                       vertices[4],vertices[7]);// Lat plane
1420  G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1421                                       vertices[7],vertices[6]);// Lat plane
1422  G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1423                                       vertices[5],vertices[6]);// Lat plane
1424  G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1425                                       vertices[6],vertices[7]);// fDz plane
1426  rand = G4UniformRand();
1427  area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
1428  chose = rand*area;
1429
1430  if ( ( chose < Surface0)
1431    || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
1432  {                                        // fDz or -fDz Plane
1433    ipl = G4int(G4UniformRand()*4);
1434    j = (ipl+1)%4;
1435    if(chose < Surface0)
1436    { 
1437      zp = -fDz;
1438      u = fVertices[ipl]; v = fVertices[j];
1439      w = fVertices[(ipl+3)%4]; 
1440    }
1441    else
1442    {
1443      zp = fDz;
1444      u = fVertices[ipl+4]; v = fVertices[j+4];
1445      w = fVertices[(ipl+3)%4+4]; 
1446    }
1447    alfa = G4UniformRand();
1448    beta = G4UniformRand();
1449    lambda1=alfa*beta;
1450    lambda0=alfa-lambda1;
1451    v = v-u;
1452    w = w-u;
1453    v = u+lambda0*v+lambda1*w;
1454  }
1455  else                                     // Lateral Plane Twisted or Not
1456  {
1457    if (chose < Surface0+Surface1) { ipl=0; }
1458    else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
1459    else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
1460    else { ipl=3; }
1461    j = (ipl+1)%4;
1462    zp = -fDz+G4UniformRand()*2*fDz;
1463    cf = 0.5*(fDz-zp)/fDz;
1464    u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
1465    v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]); 
1466    v = u+(v-u)*G4UniformRand();
1467  }
1468  point=G4ThreeVector(v.x(),v.y(),zp);
1469
1470  return point;
1471}
1472
1473// --------------------------------------------------------------------
1474
1475G4bool G4GenericTrap::ComputeIsTwisted() 
1476{
1477  // Computes tangents of twist angles (angles between projections on XY plane
1478  // of corresponding -dz +dz edges).
1479
1480  G4bool twisted = false;
1481  G4double dx1, dy1, dx2, dy2;
1482  G4int nv = fgkNofVertices/2;
1483
1484  for ( G4int i=0; i<4; i++ )
1485  {
1486    dx1 = fVertices[(i+1)%nv].x()-fVertices[i].x();
1487    dy1 = fVertices[(i+1)%nv].y()-fVertices[i].y();
1488    if ( (dx1 == 0) && (dy1 == 0) ) { continue; }
1489
1490    dx2 = fVertices[nv+(i+1)%nv].x()-fVertices[nv+i].x();
1491    dy2 = fVertices[nv+(i+1)%nv].y()-fVertices[nv+i].y();
1492
1493    if ( dx2 == 0 && dy2 == 0 ) { continue; }
1494    G4double twist_angle = std::fabs(dy1*dx2 - dx1*dy2);       
1495    if ( twist_angle < fgkTolerance ) { continue; }
1496    twisted = true;
1497    SetTwistAngle(i,twist_angle);
1498  }
1499
1500  return twisted;
1501}
1502
1503// --------------------------------------------------------------------
1504
1505G4bool G4GenericTrap::CheckOrder(const std::vector<G4TwoVector>& vertices) const
1506{
1507  // Test if the vertices are in a clockwise order, if not reorder them.
1508  // Also test if they're well defined without crossing opposite segments
1509
1510  G4bool clockwise_order=true;
1511  G4double sum1 = 0.;
1512  G4double sum2 = 0.;
1513  G4int j;
1514
1515  for (G4int i=0; i<4; i++)
1516  {
1517    j = (i+1)%4;
1518    sum1 += vertices[i].x()*vertices[j].y() - vertices[j].x()*vertices[i].y();
1519    sum2 += vertices[i+4].x()*vertices[j+4].y()
1520          - vertices[j+4].x()*vertices[i+4].y();
1521  }
1522  if (sum1*sum2 < -fgkTolerance)
1523  {
1524    G4String errorDescription = "InvalidSetup in \"";
1525    errorDescription += GetName();
1526    errorDescription += "\"";
1527
1528    G4Exception("G4GenericTrap::CheckOrder()", errorDescription, FatalException,
1529                "Lower/upper faces defined with opposite clockwise.");
1530   }
1531   
1532   if ((sum1 > 0.)||(sum2 > 0.))
1533   {
1534     G4String errorDescription = "WarningSetup in \"";
1535     errorDescription += GetName();
1536     errorDescription += "\"";
1537     G4Exception("G4GenericTrap::CheckOrder()", errorDescription, JustWarning,
1538       "Vertices must be defined in clockwise in XY planes! Re-ordering.. ");
1539     clockwise_order = false;
1540   }
1541
1542   // Check for illegal crossings
1543   //
1544   G4bool illegal_cross = false;
1545   illegal_cross = IsSegCrossing(vertices[0],vertices[1],
1546                                 vertices[2],vertices[3]);
1547   if (!illegal_cross)
1548   {
1549     illegal_cross = IsSegCrossing(vertices[4],vertices[5],
1550                                   vertices[6],vertices[7]);
1551   }
1552   if (illegal_cross)
1553   {
1554      G4String errorDescription = "InvalidSetup in \"";
1555      errorDescription += GetName();
1556      errorDescription += "\"";
1557      G4Exception("G4GenericTrap::CheckOrderAndSetup()",
1558                  errorDescription, FatalException,
1559                  "Malformed polygone with opposite sides.");
1560   }
1561   return clockwise_order;
1562}
1563
1564// --------------------------------------------------------------------
1565
1566void G4GenericTrap::ReorderVertices(std::vector<G4ThreeVector>& vertices) const
1567{
1568  // Reorder the vector of vertices
1569
1570  std::vector<G4ThreeVector> oldVertices(vertices);
1571
1572  for ( G4int i=0; i < G4int(oldVertices.size()); ++i )
1573  {
1574    vertices[i] = oldVertices[oldVertices.size()-1-i];
1575  } 
1576} 
1577 
1578// --------------------------------------------------------------------
1579
1580G4bool
1581G4GenericTrap::IsSegCrossing(const G4TwoVector& a, const G4TwoVector& b, 
1582                             const G4TwoVector& c, const G4TwoVector& d) const
1583{ 
1584  // Check if segments [A,B] and [C,D] are crossing
1585
1586  G4bool stand1 = false;
1587  G4bool stand2 = false;
1588  G4double dx1,dx2,xm=0.,ym=0.,a1=0.,a2=0.,b1=0.,b2=0.;
1589  dx1=(b-a).x();
1590  dx2=(d-c).x();
1591
1592  if( std::fabs(dx1) < fgkTolerance )  { stand1 = true; }
1593  if( std::fabs(dx2) < fgkTolerance )  { stand2 = true; }
1594  if (!stand1)
1595  {
1596    a1 = (b.x()*a.y()-a.x()*b.y())/dx1;
1597    b1 = (b-a).y()/dx1;
1598  }
1599  if (!stand2)
1600  {
1601    a2 = (d.x()*c.y()-c.x()*d.y())/dx2;
1602    b2 = (d-c).y()/dx2;
1603  }   
1604  if (stand1 && stand2)
1605  {
1606    // Segments parallel and vertical
1607    //
1608    if (std::fabs(a.x()-c.x())<fgkTolerance)
1609    {
1610       // Check if segments are overlapping
1611       //
1612       if ( ((c.y()-a.y())*(c.y()-b.y())<-fgkTolerance)
1613         || ((d.y()-a.y())*(d.y()-b.y())<-fgkTolerance)
1614         || ((a.y()-c.y())*(a.y()-d.y())<-fgkTolerance)
1615         || ((b.y()-c.y())*(b.y()-d.y())<-fgkTolerance) )  { return true; }
1616
1617       return false;
1618    }
1619    // Different x values
1620    //
1621    return false;
1622  }
1623   
1624  if (stand1)    // First segment vertical
1625  {
1626    xm = a.x();
1627    ym = a2+b2*xm; 
1628  }
1629  else
1630  {
1631    if (stand2)  // Second segment vertical
1632    {
1633      xm = c.x();
1634      ym = a1+b1*xm;
1635    }
1636    else  // Normal crossing
1637    {
1638      if (std::fabs(b1-b2) < fgkTolerance)
1639      {
1640        // Parallel segments, are they aligned
1641        //
1642        if (std::fabs(c.y()-(a1+b1*c.x())) > fgkTolerance)  { return false; }
1643
1644        // Aligned segments, are they overlapping
1645        //
1646        if ( ((c.x()-a.x())*(c.x()-b.x())<-fgkTolerance)
1647          || ((d.x()-a.x())*(d.x()-b.x())<-fgkTolerance)
1648          || ((a.x()-c.x())*(a.x()-d.x())<-fgkTolerance)
1649          || ((b.x()-c.x())*(b.x()-d.x())<-fgkTolerance) )  { return true; }
1650
1651        return false;
1652      }
1653      xm = (a1-a2)/(b2-b1);
1654      ym = (a1*b2-a2*b1)/(b2-b1);
1655    }
1656  }
1657
1658  // Check if crossing point is both between A,B and C,D
1659  //
1660  G4double check = (xm-a.x())*(xm-b.x())+(ym-a.y())*(ym-b.y());
1661  if (check > -fgkTolerance)  { return false; }
1662  check = (xm-c.x())*(xm-d.x())+(ym-c.y())*(ym-d.y());
1663  if (check > -fgkTolerance)  { return false; }
1664
1665  return true;
1666}
1667
1668// --------------------------------------------------------------------
1669
1670G4VFacet*
1671G4GenericTrap::MakeDownFacet(const std::vector<G4ThreeVector>& fromVertices, 
1672                             G4int ind1, G4int ind2, G4int ind3) const 
1673{
1674  // Create a triangular facet from the polygon points given by indices
1675  // forming the down side ( the normal goes in -z)
1676  // Do not create facet if 2 vertices are the same
1677
1678  if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1679       (fromVertices[ind2] == fromVertices[ind3]) ||
1680       (fromVertices[ind1] == fromVertices[ind3]) )  { return 0; }
1681
1682  std::vector<G4ThreeVector> vertices;
1683  vertices.push_back(fromVertices[ind1]);
1684  vertices.push_back(fromVertices[ind2]);
1685  vertices.push_back(fromVertices[ind3]);
1686 
1687  // first vertex most left
1688  //
1689  G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1690
1691  if ( cross.z() > 0.0 )
1692  {
1693    // Should not happen, as vertices should have been reordered at this stage
1694
1695    G4String errorDescription = "InvalidSetup in \"";
1696    errorDescription += GetName();
1697    errorDescription += "\"";
1698    G4Exception("G4GenericTrap::MakeDownFacet", errorDescription,
1699                FatalException, "Vertices in wrong order.");
1700  }
1701 
1702  return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
1703}
1704
1705// --------------------------------------------------------------------
1706
1707G4VFacet*
1708G4GenericTrap::MakeUpFacet(const std::vector<G4ThreeVector>& fromVertices, 
1709                           G4int ind1, G4int ind2, G4int ind3) const     
1710{
1711  // Create a triangular facet from the polygon points given by indices
1712  // forming the upper side ( z>0 )
1713
1714  // Do not create facet if 2 vertices are the same
1715  //
1716  if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1717       (fromVertices[ind2] == fromVertices[ind3]) ||
1718       (fromVertices[ind1] == fromVertices[ind3]) )  { return 0; }
1719
1720  std::vector<G4ThreeVector> vertices;
1721  vertices.push_back(fromVertices[ind1]);
1722  vertices.push_back(fromVertices[ind2]);
1723  vertices.push_back(fromVertices[ind3]);
1724 
1725  // First vertex most left
1726  //
1727  G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1728
1729  if ( cross.z() < 0.0 )
1730  {
1731    // Should not happen, as vertices should have been reordered at this stage
1732
1733    G4String errorDescription = "InvalidSetup in \"";
1734    errorDescription += GetName();
1735    errorDescription += "\"";
1736    G4Exception("G4GenericTrap::MakeUpFacet", errorDescription,
1737                FatalException, "Vertices in wrong order.");
1738  }
1739 
1740  return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
1741}     
1742
1743// --------------------------------------------------------------------
1744
1745G4VFacet*
1746G4GenericTrap::MakeSideFacet(const G4ThreeVector& downVertex0,
1747                             const G4ThreeVector& downVertex1,
1748                             const G4ThreeVector& upVertex1,
1749                             const G4ThreeVector& upVertex0) const     
1750{
1751  // Creates a triangular facet from the polygon points given by indices
1752  // forming the upper side ( z>0 )
1753
1754  if ( (downVertex0 == downVertex1) && (upVertex0 == upVertex1) )
1755  {
1756    return 0;
1757  }
1758
1759  if ( downVertex0 == downVertex1 )
1760  {
1761    return new G4TriangularFacet(downVertex0, upVertex1, upVertex0, ABSOLUTE);
1762  }
1763
1764  if ( upVertex0 == upVertex1 )
1765  {
1766    return new G4TriangularFacet(downVertex0, downVertex1, upVertex0, ABSOLUTE);
1767  }
1768
1769  return new G4QuadrangularFacet(downVertex0, downVertex1, 
1770                                 upVertex1, upVertex0, ABSOLUTE);   
1771}   
1772
1773// --------------------------------------------------------------------
1774
1775G4TessellatedSolid* G4GenericTrap::CreateTessellatedSolid() const
1776{
1777  // 3D vertices
1778  //
1779  G4int nv = fgkNofVertices/2;
1780  std::vector<G4ThreeVector> downVertices;
1781  for ( G4int i=0; i<nv; i++ )
1782  { 
1783    downVertices.push_back(G4ThreeVector(fVertices[i].x(),
1784                                         fVertices[i].y(), -fDz));
1785  }
1786
1787  std::vector<G4ThreeVector> upVertices;
1788  for ( G4int i=nv; i<2*nv; i++ )
1789  { 
1790    upVertices.push_back(G4ThreeVector(fVertices[i].x(),
1791                                       fVertices[i].y(), fDz));
1792  }
1793                                         
1794  // Reorder vertices if they are not ordered anti-clock wise
1795  //
1796  G4ThreeVector cross
1797    = (downVertices[1]-downVertices[0]).cross(downVertices[2]-downVertices[1]);
1798   G4ThreeVector cross1
1799    = (upVertices[1]-upVertices[0]).cross(upVertices[2]-upVertices[1]);
1800  if ( (cross.z() > 0.0) || (cross1.z() > 0.0) )
1801  {
1802    ReorderVertices(downVertices);
1803    ReorderVertices(upVertices);
1804  }
1805   
1806  G4TessellatedSolid* tessellatedSolid = new G4TessellatedSolid(GetName());
1807 
1808  G4VFacet* facet = 0;
1809  facet = MakeDownFacet(downVertices, 0, 1, 2);
1810  if (facet)  { tessellatedSolid->AddFacet( facet ); }
1811  facet = MakeDownFacet(downVertices, 0, 2, 3);
1812  if (facet)  { tessellatedSolid->AddFacet( facet ); }
1813  facet = MakeUpFacet(upVertices, 0, 2, 1);
1814  if (facet)  { tessellatedSolid->AddFacet( facet ); }
1815  facet = MakeUpFacet(upVertices, 0, 3, 2);
1816  if (facet)  { tessellatedSolid->AddFacet( facet ); }
1817
1818  // The quadrangular sides
1819  //
1820  for ( G4int i = 0; i < nv; ++i )
1821  {
1822    G4int j = (i+1) % nv;
1823    facet = MakeSideFacet(downVertices[j], downVertices[i], 
1824                          upVertices[i], upVertices[j]);
1825
1826    if ( facet )  { tessellatedSolid->AddFacet( facet ); }
1827  }
1828
1829  tessellatedSolid->SetSolidClosed(true);
1830
1831  return tessellatedSolid;
1832} 
1833
1834// --------------------------------------------------------------------
1835
1836void G4GenericTrap::ComputeBBox() 
1837{
1838  // Computes bounding box for a shape.
1839
1840  G4double minX, maxX, minY, maxY;
1841  minX = maxX = fVertices[0].x();
1842  minY = maxY = fVertices[0].y();
1843   
1844  for (G4int i=1; i< fgkNofVertices; i++)
1845  {
1846    if (minX>fVertices[i].x()) { minX=fVertices[i].x(); }
1847    if (maxX<fVertices[i].x()) { maxX=fVertices[i].x(); }
1848    if (minY>fVertices[i].y()) { minY=fVertices[i].y(); }
1849    if (maxY<fVertices[i].y()) { maxY=fVertices[i].y(); }
1850  }
1851  fMinBBoxVector = G4ThreeVector(minX,minY,-fDz);
1852  fMaxBBoxVector = G4ThreeVector(maxX,maxY, fDz);
1853}
1854
1855// --------------------------------------------------------------------
1856
1857G4Polyhedron* G4GenericTrap::GetPolyhedron () const
1858{
1859
1860#ifdef G4TESS_TEST
1861  if ( fTessellatedSolid )
1862  { 
1863    return fTessellatedSolid->GetPolyhedron();
1864  }
1865#endif 
1866 
1867  if ( (!fpPolyhedron)
1868    || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1869        fpPolyhedron->GetNumberOfRotationSteps()) )
1870  {
1871    delete fpPolyhedron;
1872    fpPolyhedron = CreatePolyhedron();
1873  }
1874  return fpPolyhedron;
1875}   
1876
1877// --------------------------------------------------------------------
1878
1879void G4GenericTrap::DescribeYourselfTo(G4VGraphicsScene& scene) const
1880{
1881
1882#ifdef G4TESS_TEST
1883  if ( fTessellatedSolid )
1884  { 
1885    return fTessellatedSolid->DescribeYourselfTo(scene);
1886  }
1887#endif 
1888 
1889  scene.AddSolid(*this);
1890}
1891
1892// --------------------------------------------------------------------
1893
1894G4VisExtent G4GenericTrap::GetExtent() const
1895{
1896  // Computes bounding vectors for the shape
1897
1898#ifdef G4TESS_TEST
1899  if ( fTessellatedSolid )
1900  { 
1901    return fTessellatedSolid->GetExtent();
1902  } 
1903#endif
1904   
1905  G4double Dx,Dy;
1906  G4ThreeVector minVec = GetMinimumBBox();
1907  G4ThreeVector maxVec = GetMaximumBBox();
1908  Dx = 0.5*(maxVec.x()- minVec.y());
1909  Dy = 0.5*(maxVec.y()- minVec.y());
1910
1911  return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz); 
1912}   
1913
1914// --------------------------------------------------------------------
1915
1916G4Polyhedron* G4GenericTrap::CreatePolyhedron() const
1917{
1918
1919#ifdef G4TESS_TEST
1920  if ( fTessellatedSolid )
1921  { 
1922    return fTessellatedSolid->CreatePolyhedron();
1923  } 
1924#endif
1925 
1926  // Approximation of Twisted Side
1927  // Construct extra Points, if Twisted Side
1928  //
1929  G4PolyhedronArbitrary* polyhedron;
1930  size_t nVertices, nFacets;
1931
1932  G4int subdivisions=0;
1933  G4int i;
1934  if(fIsTwisted)
1935  {
1936    if ( GetVisSubdivisions()!= 0 )
1937    {
1938      subdivisions=GetVisSubdivisions();
1939    }
1940    else
1941    {
1942      // Estimation of Number of Subdivisions for smooth visualisation
1943      //
1944      G4double maxTwist=0.;
1945      for(i=0; i<4; i++)
1946      {
1947        if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
1948      }
1949
1950      // Computes bounding vectors for the shape
1951      //
1952      G4double Dx,Dy;
1953      G4ThreeVector minVec = GetMinimumBBox();
1954      G4ThreeVector maxVec = GetMaximumBBox();
1955      Dx = 0.5*(maxVec.x()- minVec.y());
1956      Dy = 0.5*(maxVec.y()- minVec.y());
1957      if (Dy > Dx)  { Dx=Dy; }
1958   
1959      subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
1960      if (subdivisions<4)  { subdivisions=4; }
1961      if (subdivisions>30) { subdivisions=30; }
1962    }
1963  }
1964  G4int sub4=4*subdivisions;
1965  nVertices = 8+subdivisions*4;
1966  nFacets = 6+subdivisions*4;
1967  G4double cf=1./(subdivisions+1);
1968  polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
1969
1970  // Add Vertex
1971  //
1972  for (i=0;i<4;i++)
1973  {
1974    polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
1975                                        fVertices[i].y(),-fDz));
1976  }
1977  for( i=0;i<subdivisions;i++)
1978  {
1979    for(G4int j=0;j<4;j++)
1980    {
1981      G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
1982      polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
1983    }   
1984  }
1985  for (i=4;i<8;i++)
1986  {
1987    polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
1988                                        fVertices[i].y(),fDz));
1989  }
1990
1991  // Add Facets
1992  //
1993  polyhedron->AddFacet(1,4,3,2);  //Z-plane
1994  for (i=0;i<subdivisions+1;i++)
1995  {
1996    G4int is=i*4;
1997    polyhedron->AddFacet(5+is,8+is,4+is,1+is);
1998    polyhedron->AddFacet(8+is,7+is,3+is,4+is);
1999    polyhedron->AddFacet(7+is,6+is,2+is,3+is);
2000    polyhedron->AddFacet(6+is,5+is,1+is,2+is); 
2001  }
2002  polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4);  //Z-plane
2003
2004  return (G4Polyhedron*) polyhedron;
2005}
2006
2007// --------------------------------------------------------------------
2008
2009G4NURBS*  G4GenericTrap::CreateNURBS() const
2010{
2011#ifdef G4TESS_TEST
2012  if ( fTessellatedSolid )
2013  { 
2014    return fTessellatedSolid->CreateNURBS();
2015  }
2016#endif
2017
2018  // Computes bounding vectors for the shape
2019  //
2020  G4double Dx,Dy;
2021  G4ThreeVector minVec = GetMinimumBBox();
2022  G4ThreeVector maxVec = GetMaximumBBox();
2023  Dx = 0.5*(maxVec.x()- minVec.y());
2024  Dy = 0.5*(maxVec.y()- minVec.y());
2025
2026  return new G4NURBSbox (Dx, Dy, fDz);
2027}   
2028
2029// --------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.