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