Geant4  10.02
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 83851 2014-09-19 10:12:12Z 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  G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
509  JustWarning, "Point p is not on surface !?" );
510  sumnorm=apprnorm;
511  // Add Approximative Surface Normal Calculation?
512  }
513  else if ( noSurfaces == 1 ) { sumnorm = sumnorm; }
514  else { sumnorm = sumnorm.unit(); }
515 
516  return sumnorm ;
517 }
518 
519 // --------------------------------------------------------------------
520 
522  const G4int ipl ) const
523 {
524  // Return normal to given lateral plane ipl
525 
526 #ifdef G4TESS_TEST
527  if ( fTessellatedSolid )
528  {
529  return fTessellatedSolid->SurfaceNormal(p);
530  }
531 #endif
532 
533  G4ThreeVector lnorm, norm(0.,0.,0.), p0,p1,p2;
534 
535  G4double distz = fDz-p.z();
536  G4int i=ipl; // current plane index
537 
538  G4TwoVector u,v;
539  G4ThreeVector r1,r2,r3,r4;
540  G4double cf = 0.5*(fDz-p.z())/fDz;
541  G4int j=(i+1)%4;
542 
543  u=fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]);
544  v=fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
545 
546  // Compute cross product
547  //
548  p0=G4ThreeVector(u.x(),u.y(),p.z());
549 
550  if (std::fabs(distz)<halfCarTolerance)
551  {
552  p1=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);distz=-1;}
553  else
554  {
555  p1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
556  }
557  p2=G4ThreeVector(v.x(),v.y(),p.z());
558 
559  // Collapsed vertices
560  //
561  if ( (p2-p0).mag2() < kCarTolerance )
562  {
563  if ( std::fabs(p.z()+fDz) > halfCarTolerance )
564  {
565  p2=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
566  }
567  else
568  {
569  p2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
570  }
571  }
572  lnorm=-(p1-p0).cross(p2-p0);
573  if (distz>-halfCarTolerance) { lnorm=-lnorm.unit(); }
574  else { lnorm=lnorm.unit(); }
575 
576  // Adjust Normal for Twisted Surface
577  //
578  if( (fIsTwisted) && (GetTwistAngle(ipl)!=0) )
579  {
580  G4double normP=(p2-p0).mag();
581  if(normP)
582  {
583  G4double proj=(p-p0).dot(p2-p0)/normP;
584  if (proj<0) { proj=0; }
585  if (proj>normP) { proj=normP; }
586 
587  r1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
588  r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
589  r3=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);
590  r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
591  r1=r1+proj*(r2-r1)/normP;
592  r3=r3+proj*(r4-r3)/normP;
593  r2=r1-r3;
594  r4=r2.cross(p2-p0);r4=r4.unit();
595  lnorm=r4;
596  }
597  } // End if fIsTwisted
598 
599  return lnorm;
600 }
601 
602 // --------------------------------------------------------------------
603 
605  const G4ThreeVector& v,
606  const G4int ipl) const
607 {
608  // Computes distance to plane ipl :
609  // ipl=0 : points 0,4,1,5
610  // ipl=1 : points 1,5,2,6
611  // ipl=2 : points 2,6,3,7
612  // ipl=3 : points 3,7,0,4
613 
614  G4double xa,xb,xc,xd,ya,yb,yc,yd;
615 
616  G4int j = (ipl+1)%4;
617 
618  xa=fVertices[ipl].x();
619  ya=fVertices[ipl].y();
620  xb=fVertices[ipl+4].x();
621  yb=fVertices[ipl+4].y();
622  xc=fVertices[j].x();
623  yc=fVertices[j].y();
624  xd=fVertices[4+j].x();
625  yd=fVertices[4+j].y();
626 
627  G4double dz2 =0.5/fDz;
628  G4double tx1 =dz2*(xb-xa);
629  G4double ty1 =dz2*(yb-ya);
630  G4double tx2 =dz2*(xd-xc);
631  G4double ty2 =dz2*(yd-yc);
632  G4double dzp =fDz+p.z();
633  G4double xs1 =xa+tx1*dzp;
634  G4double ys1 =ya+ty1*dzp;
635  G4double xs2 =xc+tx2*dzp;
636  G4double ys2 =yc+ty2*dzp;
637  G4double dxs =xs2-xs1;
638  G4double dys =ys2-ys1;
639  G4double dtx =tx2-tx1;
640  G4double dty =ty2-ty1;
641 
642  G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
643  G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
644  + tx1*ys2-tx2*ys1)*v.z();
645  G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
647  G4double x1,x2,y1,y2,xp,yp,zi;
648 
649  if (std::fabs(a)<kCarTolerance)
650  {
651  if (std::fabs(b)<kCarTolerance) { return kInfinity; }
652  q=-c/b;
653 
654  // Check if Point is on the Surface
655 
656  if (q>-halfCarTolerance)
657  {
658  if (q<halfCarTolerance)
659  {
660  if (NormalToPlane(p,ipl).dot(v)<=0)
661  { if(Inside(p) != kOutside) { return 0.; } }
662  else
663  { return kInfinity; }
664  }
665 
666  // Check the Intersection
667  //
668  zi=p.z()+q*v.z();
669  if (std::fabs(zi)<fDz)
670  {
671  x1=xs1+tx1*v.z()*q;
672  x2=xs2+tx2*v.z()*q;
673  xp=p.x()+q*v.x();
674  y1=ys1+ty1*v.z()*q;
675  y2=ys2+ty2*v.z()*q;
676  yp=p.y()+q*v.y();
677  zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
678  if (zi<=halfCarTolerance) { return q; }
679  }
680  }
681  return kInfinity;
682  }
683  G4double d=b*b-4*a*c;
684  if (d>=0)
685  {
686  if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
687  else { q=0.5*(-b+std::sqrt(d))/a; }
688 
689  // Check if Point is on the Surface
690  //
691  if (q>-halfCarTolerance)
692  {
693  if(q<halfCarTolerance)
694  {
695  if (NormalToPlane(p,ipl).dot(v)<=0)
696  {
697  if(Inside(p)!= kOutside) { return 0.; }
698  }
699  else // Check second root; return kInfinity
700  {
701  if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
702  else { q=0.5*(-b-std::sqrt(d))/a; }
703  if (q<=halfCarTolerance) { return kInfinity; }
704  }
705  }
706  // Check the Intersection
707  //
708  zi=p.z()+q*v.z();
709  if (std::fabs(zi)<fDz)
710  {
711  x1=xs1+tx1*v.z()*q;
712  x2=xs2+tx2*v.z()*q;
713  xp=p.x()+q*v.x();
714  y1=ys1+ty1*v.z()*q;
715  y2=ys2+ty2*v.z()*q;
716  yp=p.y()+q*v.y();
717  zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
718  if (zi<=halfCarTolerance) { return q; }
719  }
720  }
721  if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
722  else { q=0.5*(-b-std::sqrt(d))/a; }
723 
724  // Check if Point is on the Surface
725  //
726  if (q>-halfCarTolerance)
727  {
728  if(q<halfCarTolerance)
729  {
730  if (NormalToPlane(p,ipl).dot(v)<=0)
731  {
732  if(Inside(p) != kOutside) { return 0.; }
733  }
734  else // Check second root; return kInfinity.
735  {
736  if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
737  else { q=0.5*(-b+std::sqrt(d))/a; }
738  if (q<=halfCarTolerance) { return kInfinity; }
739  }
740  }
741  // Check the Intersection
742  //
743  zi=p.z()+q*v.z();
744  if (std::fabs(zi)<fDz)
745  {
746  x1=xs1+tx1*v.z()*q;
747  x2=xs2+tx2*v.z()*q;
748  xp=p.x()+q*v.x();
749  y1=ys1+ty1*v.z()*q;
750  y2=ys2+ty2*v.z()*q;
751  yp=p.y()+q*v.y();
752  zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
753  if (zi<=halfCarTolerance) { return q; }
754  }
755  }
756  }
757  return kInfinity;
758 }
759 
760 // --------------------------------------------------------------------
761 
763  const G4ThreeVector& v) const
764 {
765 #ifdef G4TESS_TEST
766  if ( fTessellatedSolid )
767  {
768  return fTessellatedSolid->DistanceToIn(p, v);
769  }
770 #endif
771 
772  G4double dist[5];
774 
775  // Check lateral faces
776  //
777  G4int i;
778  for (i=0; i<4; i++)
779  {
780  dist[i]=DistToPlane(p, v, i);
781  }
782 
783  // Check Z planes
784  //
785  dist[4]=kInfinity;
786  if (std::fabs(p.z())>fDz-halfCarTolerance)
787  {
788  if (v.z())
789  {
790  G4ThreeVector pt;
791  if (p.z()>0)
792  {
793  dist[4] = (fDz-p.z())/v.z();
794  }
795  else
796  {
797  dist[4] = (-fDz-p.z())/v.z();
798  }
799  if (dist[4]<-halfCarTolerance)
800  {
801  dist[4]=kInfinity;
802  }
803  else
804  {
805  if(dist[4]<halfCarTolerance)
806  {
807  if(p.z()>0) { n=G4ThreeVector(0,0,1); }
808  else { n=G4ThreeVector(0,0,-1); }
809  if (n.dot(v)<0) { dist[4]=0.; }
810  else { dist[4]=kInfinity; }
811  }
812  pt=p+dist[4]*v;
813  if (Inside(pt)==kOutside) { dist[4]=kInfinity; }
814  }
815  }
816  }
817  G4double distmin = dist[0];
818  for (i=1;i<5;i++)
819  {
820  if (dist[i] < distmin) { distmin = dist[i]; }
821  }
822 
823  if (distmin<halfCarTolerance) { distmin=0.; }
824 
825  return distmin;
826 }
827 
828 // --------------------------------------------------------------------
829 
831 {
832  // Computes the closest distance from given point to this shape
833 
834 #ifdef G4TESS_TEST
835  if ( fTessellatedSolid )
836  {
837  return fTessellatedSolid->DistanceToIn(p);
838  }
839 #endif
840 
841  G4double safz = std::fabs(p.z())-fDz;
842  if(safz<0) { safz=0; }
843 
844  G4int iseg;
845  G4double safe = safz;
846  G4double safxy = safz;
847 
848  for (iseg=0; iseg<4; iseg++)
849  {
850  safxy = SafetyToFace(p,iseg);
851  if (safxy>safe) { safe=safxy; }
852  }
853 
854  return safe;
855 }
856 
857 // --------------------------------------------------------------------
858 
859 G4double
861 {
862  // Estimate distance to lateral plane defined by segment iseg in range [0,3]
863  // Might be negative: plane seen only from inside
864 
865  G4ThreeVector p1,norm;
866  G4double safe;
867 
868  p1=G4ThreeVector(fVertices[iseg].x(),fVertices[iseg].y(),-fDz);
869 
870  norm=NormalToPlane(p,iseg);
871  safe = (p-p1).dot(norm); // Can be negative
872 
873  return safe;
874 }
875 
876 // --------------------------------------------------------------------
877 
878 G4double
880  const G4ThreeVector& v, const G4int ipl) const
881 {
882  G4double xa=fVertices[ipl].x();
883  G4double ya=fVertices[ipl].y();
884  G4double xb=fVertices[ipl+4].x();
885  G4double yb=fVertices[ipl+4].y();
886  G4int j=(ipl+1)%4;
887  G4double xc=fVertices[j].x();
888  G4double yc=fVertices[j].y();
889  G4double zab=2*fDz;
890  G4double zac=0;
891 
892  if ( (std::fabs(xa-xc)+std::fabs(ya-yc)) < halfCarTolerance )
893  {
894  xc=fVertices[j+4].x();
895  yc=fVertices[j+4].y();
896  zac=2*fDz;
897  zab=2*fDz;
898 
899  //Line case
900  //
901  if ( (std::fabs(xb-xc)+std::fabs(yb-yc)) < halfCarTolerance )
902  {
903  return kInfinity;
904  }
905  }
906  G4double a=(yb-ya)*zac-(yc-ya)*zab;
907  G4double b=(xc-xa)*zab-(xb-xa)*zac;
908  G4double c=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
909  G4double d=-xa*a-ya*b+fDz*c;
910  G4double t=a*v.x()+b*v.y()+c*v.z();
911 
912  if (t!=0)
913  {
914  t=-(a*p.x()+b*p.y()+c*p.z()+d)/t;
915  }
916  if ( (t<halfCarTolerance) && (t>-halfCarTolerance) )
917  {
918  if (NormalToPlane(p,ipl).dot(v)<kCarTolerance)
919  {
920  t=kInfinity;
921  }
922  else
923  {
924  t=0;
925  }
926  }
927  if (Inside(p+v*t) != kSurface) { t=kInfinity; }
928 
929  return t;
930 }
931 
932 // --------------------------------------------------------------------
933 
935  const G4ThreeVector& v,
936  const G4bool calcNorm,
937  G4bool* validNorm,
938  G4ThreeVector* n) const
939 {
940 #ifdef G4TESS_TEST
941  if ( fTessellatedSolid )
942  {
943  return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
944  }
945 #endif
946 
947  G4double distmin;
948  G4bool lateral_cross = false;
949  ESide side = kUndefined;
950 
951  if (calcNorm) { *validNorm=true; } // All normals are valid
952 
953  if (v.z() < 0)
954  {
955  distmin=(-fDz-p.z())/v.z();
956  if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
957  }
958  else
959  {
960  if (v.z() > 0)
961  {
962  distmin = (fDz-p.z())/v.z();
963  if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); }
964  }
965  else { distmin = kInfinity; }
966  }
967 
968  G4double dz2 =0.5/fDz;
969  G4double xa,xb,xc,xd;
970  G4double ya,yb,yc,yd;
971 
972  for (G4int ipl=0; ipl<4; ipl++)
973  {
974  G4int j = (ipl+1)%4;
975  xa=fVertices[ipl].x();
976  ya=fVertices[ipl].y();
977  xb=fVertices[ipl+4].x();
978  yb=fVertices[ipl+4].y();
979  xc=fVertices[j].x();
980  yc=fVertices[j].y();
981  xd=fVertices[4+j].x();
982  yd=fVertices[4+j].y();
983 
984  if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
985  || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
986  {
987  G4double q=DistToTriangle(p,v,ipl) ;
988  if ( (q>=0) && (q<distmin) )
989  {
990  distmin=q;
991  lateral_cross=true;
992  side=ESide(ipl+1);
993  }
994  continue;
995  }
996  G4double tx1 =dz2*(xb-xa);
997  G4double ty1 =dz2*(yb-ya);
998  G4double tx2 =dz2*(xd-xc);
999  G4double ty2 =dz2*(yd-yc);
1000  G4double dzp =fDz+p.z();
1001  G4double xs1 =xa+tx1*dzp;
1002  G4double ys1 =ya+ty1*dzp;
1003  G4double xs2 =xc+tx2*dzp;
1004  G4double ys2 =yc+ty2*dzp;
1005  G4double dxs =xs2-xs1;
1006  G4double dys =ys2-ys1;
1007  G4double dtx =tx2-tx1;
1008  G4double dty =ty2-ty1;
1009  G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
1010  G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
1011  + tx1*ys2-tx2*ys1)*v.z();
1012  G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
1013  G4double q=kInfinity;
1014 
1015  if (std::fabs(a) < kCarTolerance)
1016  {
1017  if (std::fabs(b) < kCarTolerance) { continue; }
1018  q=-c/b;
1019 
1020  // Check for Point on the Surface
1021  //
1022  if ((q > -halfCarTolerance) && (q < distmin))
1023  {
1024  if (q < halfCarTolerance)
1025  {
1026  if (NormalToPlane(p,ipl).dot(v)<0.) { continue; }
1027  }
1028  distmin =q;
1029  lateral_cross=true;
1030  side=ESide(ipl+1);
1031  }
1032  continue;
1033  }
1034  G4double d=b*b-4*a*c;
1035  if (d >= 0.)
1036  {
1037  if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1038  else { q=0.5*(-b+std::sqrt(d))/a; }
1039 
1040  // Check for Point on the Surface
1041  //
1042  if (q > -halfCarTolerance )
1043  {
1044  if (q < distmin)
1045  {
1046  if(q < halfCarTolerance)
1047  {
1048  if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
1049  {
1050  if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1051  else { q=0.5*(-b-std::sqrt(d))/a; }
1052  if (( q > halfCarTolerance) && (q < distmin))
1053  {
1054  distmin=q;
1055  lateral_cross = true;
1056  side=ESide(ipl+1);
1057  }
1058  continue;
1059  }
1060  }
1061  distmin = q;
1062  lateral_cross = true;
1063  side=ESide(ipl+1);
1064  }
1065  }
1066  else
1067  {
1068  if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1069  else { q=0.5*(-b-std::sqrt(d))/a; }
1070 
1071  // Check for Point on the Surface
1072  //
1073  if ((q > -halfCarTolerance) && (q < distmin))
1074  {
1075  if (q < halfCarTolerance)
1076  {
1077  if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
1078  {
1079  if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1080  else { q=0.5*(-b+std::sqrt(d))/a; }
1081  if ( ( q > halfCarTolerance) && (q < distmin) )
1082  {
1083  distmin=q;
1084  lateral_cross = true;
1085  side=ESide(ipl+1);
1086  }
1087  continue;
1088  }
1089  }
1090  distmin =q;
1091  lateral_cross = true;
1092  side=ESide(ipl+1);
1093  }
1094  }
1095  }
1096  }
1097  if (!lateral_cross) // Make sure that track crosses the top or bottom
1098  {
1099  if (distmin >= kInfinity) { distmin=kCarTolerance; }
1100  G4ThreeVector pt=p+distmin*v;
1101 
1102  // Check if propagated point is in the polygon
1103  //
1104  G4int i=0;
1105  if (v.z()>0.) { i=4; }
1106  std::vector<G4TwoVector> xy;
1107  for ( G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
1108 
1109  // Check Inside
1110  //
1111  if (InsidePolygone(pt,xy)==kOutside)
1112  {
1113  if(calcNorm)
1114  {
1115  if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
1116  else { side=kMZ; *n = G4ThreeVector(0,0,-1);}
1117  }
1118  return 0.;
1119  }
1120  else
1121  {
1122  if(v.z()>0) {side=kPZ;}
1123  else {side=kMZ;}
1124  }
1125  }
1126 
1127  if (calcNorm)
1128  {
1129  G4ThreeVector pt=p+v*distmin;
1130  switch (side)
1131  {
1132  case kXY0:
1133  *n=NormalToPlane(pt,0);
1134  break;
1135  case kXY1:
1136  *n=NormalToPlane(pt,1);
1137  break;
1138  case kXY2:
1139  *n=NormalToPlane(pt,2);
1140  break;
1141  case kXY3:
1142  *n=NormalToPlane(pt,3);
1143  break;
1144  case kPZ:
1145  *n=G4ThreeVector(0,0,1);
1146  break;
1147  case kMZ:
1148  *n=G4ThreeVector(0,0,-1);
1149  break;
1150  default:
1151  DumpInfo();
1152  std::ostringstream message;
1153  G4int oldprc = message.precision(16);
1154  message << "Undefined side for valid surface normal to solid." << G4endl
1155  << "Position:" << G4endl
1156  << " p.x() = " << p.x()/mm << " mm" << G4endl
1157  << " p.y() = " << p.y()/mm << " mm" << G4endl
1158  << " p.z() = " << p.z()/mm << " mm" << G4endl
1159  << "Direction:" << G4endl
1160  << " v.x() = " << v.x() << G4endl
1161  << " v.y() = " << v.y() << G4endl
1162  << " v.z() = " << v.z() << G4endl
1163  << "Proposed distance :" << G4endl
1164  << " distmin = " << distmin/mm << " mm";
1165  message.precision(oldprc);
1166  G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
1167  "GeomSolids1002", JustWarning, message);
1168  break;
1169  }
1170  }
1171 
1172  if (distmin<halfCarTolerance) { distmin=0.; }
1173 
1174  return distmin;
1175 }
1176 
1177 // --------------------------------------------------------------------
1178 
1180 {
1181 
1182 #ifdef G4TESS_TEST
1183  if ( fTessellatedSolid )
1184  {
1185  return fTessellatedSolid->DistanceToOut(p);
1186  }
1187 #endif
1188 
1189  G4double safz = fDz-std::fabs(p.z());
1190  if (safz<0) { safz = 0; }
1191 
1192  G4double safe = safz;
1193  G4double safxy = safz;
1194 
1195  for (G4int iseg=0; iseg<4; iseg++)
1196  {
1197  safxy = std::fabs(SafetyToFace(p,iseg));
1198  if (safxy < safe) { safe = safxy; }
1199  }
1200 
1201  return safe;
1202 }
1203 
1204 // --------------------------------------------------------------------
1205 
1207  const G4VoxelLimits& pVoxelLimit,
1208  const G4AffineTransform& pTransform,
1209  G4double& pMin, G4double& pMax) const
1210 {
1211 #ifdef G4TESS_TEST
1212  if ( fTessellatedSolid )
1213  {
1214  return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
1215  pTransform, pMin, pMax);
1216  }
1217 #endif
1218 
1219  // Computes bounding vectors for a shape
1220  //
1221  G4double Dx,Dy;
1222  G4ThreeVector minVec = GetMinimumBBox();
1223  G4ThreeVector maxVec = GetMaximumBBox();
1224  Dx = 0.5*(maxVec.x()- minVec.x());
1225  Dy = 0.5*(maxVec.y()- minVec.y());
1226 
1227  if (!pTransform.IsRotated())
1228  {
1229  // Special case handling for unrotated shapes
1230  // Compute x/y/z mins and maxs respecting limits, with early returns
1231  // if outside limits. Then switch() on pAxis
1232  //
1233  G4double xoffset,xMin,xMax;
1234  G4double yoffset,yMin,yMax;
1235  G4double zoffset,zMin,zMax;
1236 
1237  xoffset=pTransform.NetTranslation().x();
1238  xMin=xoffset-Dx;
1239  xMax=xoffset+Dx;
1240  if (pVoxelLimit.IsXLimited())
1241  {
1242  if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
1243  || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
1244  {
1245  return false;
1246  }
1247  else
1248  {
1249  if (xMin<pVoxelLimit.GetMinXExtent())
1250  {
1251  xMin=pVoxelLimit.GetMinXExtent();
1252  }
1253  if (xMax>pVoxelLimit.GetMaxXExtent())
1254  {
1255  xMax=pVoxelLimit.GetMaxXExtent();
1256  }
1257  }
1258  }
1259 
1260  yoffset=pTransform.NetTranslation().y();
1261  yMin=yoffset-Dy;
1262  yMax=yoffset+Dy;
1263  if (pVoxelLimit.IsYLimited())
1264  {
1265  if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
1266  || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
1267  {
1268  return false;
1269  }
1270  else
1271  {
1272  if (yMin<pVoxelLimit.GetMinYExtent())
1273  {
1274  yMin=pVoxelLimit.GetMinYExtent();
1275  }
1276  if (yMax>pVoxelLimit.GetMaxYExtent())
1277  {
1278  yMax=pVoxelLimit.GetMaxYExtent();
1279  }
1280  }
1281  }
1282 
1283  zoffset=pTransform.NetTranslation().z();
1284  zMin=zoffset-fDz;
1285  zMax=zoffset+fDz;
1286  if (pVoxelLimit.IsZLimited())
1287  {
1288  if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
1289  || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
1290  {
1291  return false;
1292  }
1293  else
1294  {
1295  if (zMin<pVoxelLimit.GetMinZExtent())
1296  {
1297  zMin=pVoxelLimit.GetMinZExtent();
1298  }
1299  if (zMax>pVoxelLimit.GetMaxZExtent())
1300  {
1301  zMax=pVoxelLimit.GetMaxZExtent();
1302  }
1303  }
1304  }
1305 
1306  switch (pAxis)
1307  {
1308  case kXAxis:
1309  pMin = xMin;
1310  pMax = xMax;
1311  break;
1312  case kYAxis:
1313  pMin = yMin;
1314  pMax = yMax;
1315  break;
1316  case kZAxis:
1317  pMin = zMin;
1318  pMax = zMax;
1319  break;
1320  default:
1321  break;
1322  }
1323  pMin-=kCarTolerance;
1324  pMax+=kCarTolerance;
1325 
1326  return true;
1327  }
1328  else
1329  {
1330  // General rotated case - create and clip mesh to boundaries
1331 
1332  G4bool existsAfterClip=false;
1333  G4ThreeVectorList *vertices;
1334 
1335  pMin=+kInfinity;
1336  pMax=-kInfinity;
1337 
1338  // Calculate rotated vertex coordinates
1339  //
1340  vertices=CreateRotatedVertices(pTransform);
1341  ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
1342  ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
1343  ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
1344 
1345  if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
1346  {
1347  existsAfterClip=true;
1348 
1349  // Add 2*tolerance to avoid precision troubles
1350  //
1351  pMin-=kCarTolerance;
1352  pMax+=kCarTolerance;
1353  }
1354  else
1355  {
1356  // Check for case where completely enveloping clipping volume.
1357  // If point inside then we are confident that the solid completely
1358  // envelopes the clipping volume. Hence set min/max extents according
1359  // to clipping volume extents along the specified axis.
1360  //
1361  G4ThreeVector clipCentre(
1362  (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
1363  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
1364  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
1365 
1366  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
1367  {
1368  existsAfterClip=true;
1369  pMin=pVoxelLimit.GetMinExtent(pAxis);
1370  pMax=pVoxelLimit.GetMaxExtent(pAxis);
1371  }
1372  }
1373  delete vertices;
1374  return existsAfterClip;
1375  }
1376 }
1377 
1378 // --------------------------------------------------------------------
1379 
1382 {
1383  // Create a List containing the transformed vertices
1384  // Ordering [0-3] -fDz cross section
1385  // [4-7] +fDz cross section such that [0] is below [4],
1386  // [1] below [5] etc.
1387  // Note: caller has deletion responsibility
1388 
1389  G4ThreeVector Min = GetMinimumBBox();
1390  G4ThreeVector Max = GetMaximumBBox();
1391 
1392  G4ThreeVectorList *vertices;
1393  vertices=new G4ThreeVectorList();
1394 
1395  if (vertices)
1396  {
1397  vertices->reserve(8);
1398  G4ThreeVector vertex0(Min.x(),Min.y(),Min.z());
1399  G4ThreeVector vertex1(Max.x(),Min.y(),Min.z());
1400  G4ThreeVector vertex2(Max.x(),Max.y(),Min.z());
1401  G4ThreeVector vertex3(Min.x(),Max.y(),Min.z());
1402  G4ThreeVector vertex4(Min.x(),Min.y(),Max.z());
1403  G4ThreeVector vertex5(Max.x(),Min.y(),Max.z());
1404  G4ThreeVector vertex6(Max.x(),Max.y(),Max.z());
1405  G4ThreeVector vertex7(Min.x(),Max.y(),Max.z());
1406 
1407  vertices->push_back(pTransform.TransformPoint(vertex0));
1408  vertices->push_back(pTransform.TransformPoint(vertex1));
1409  vertices->push_back(pTransform.TransformPoint(vertex2));
1410  vertices->push_back(pTransform.TransformPoint(vertex3));
1411  vertices->push_back(pTransform.TransformPoint(vertex4));
1412  vertices->push_back(pTransform.TransformPoint(vertex5));
1413  vertices->push_back(pTransform.TransformPoint(vertex6));
1414  vertices->push_back(pTransform.TransformPoint(vertex7));
1415  }
1416  else
1417  {
1418  G4Exception("G4GenericTrap::CreateRotatedVertices()", "FatalError",
1419  FatalException, "Out of memory - Cannot allocate vertices!");
1420  }
1421  return vertices;
1422 }
1423 
1424 // --------------------------------------------------------------------
1425 
1427 {
1428  return G4String("G4GenericTrap");
1429 }
1430 
1431 // --------------------------------------------------------------------
1432 
1434 {
1435  return new G4GenericTrap(*this);
1436 }
1437 
1438 // --------------------------------------------------------------------
1439 
1440 std::ostream& G4GenericTrap::StreamInfo(std::ostream& os) const
1441 {
1442  G4int oldprc = os.precision(16);
1443  os << "-----------------------------------------------------------\n"
1444  << " *** Dump for solid - " << GetName() << " *** \n"
1445  << " =================================================== \n"
1446  << " Solid geometry type: " << GetEntityType() << G4endl
1447  << " half length Z: " << fDz/mm << " mm \n"
1448  << " list of vertices:\n";
1449 
1450  for ( G4int i=0; i<fgkNofVertices; ++i )
1451  {
1452  os << std::setw(5) << "#" << i
1453  << " vx = " << fVertices[i].x()/mm << " mm"
1454  << " vy = " << fVertices[i].y()/mm << " mm" << G4endl;
1455  }
1456  os.precision(oldprc);
1457 
1458  return os;
1459 }
1460 
1461 // --------------------------------------------------------------------
1462 
1464 {
1465 
1466 #ifdef G4TESS_TEST
1467  if ( fTessellatedSolid )
1468  {
1470  }
1471 #endif
1472 
1473  G4ThreeVector point;
1474  G4TwoVector u,v,w;
1475  G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
1476  G4int ipl,j;
1477 
1478  std::vector<G4ThreeVector> vertices;
1479  for (G4int i=0; i<4;i++)
1480  {
1481  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1482  }
1483  for (G4int i=4; i<8;i++)
1484  {
1485  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1486  }
1487 
1488  // Surface Area of Planes(only estimation for twisted)
1489  //
1490  G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1491  vertices[2],vertices[3]);//-fDz plane
1492  G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1493  vertices[5],vertices[4]);// Lat plane
1494  G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1495  vertices[4],vertices[7]);// Lat plane
1496  G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1497  vertices[7],vertices[6]);// Lat plane
1498  G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1499  vertices[5],vertices[6]);// Lat plane
1500  G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1501  vertices[6],vertices[7]);// fDz plane
1502  rand = G4UniformRand();
1503  area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
1504  chose = rand*area;
1505 
1506  if ( ( chose < Surface0)
1507  || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
1508  { // fDz or -fDz Plane
1509  ipl = G4int(G4UniformRand()*4);
1510  j = (ipl+1)%4;
1511  if(chose < Surface0)
1512  {
1513  zp = -fDz;
1514  u = fVertices[ipl]; v = fVertices[j];
1515  w = fVertices[(ipl+3)%4];
1516  }
1517  else
1518  {
1519  zp = fDz;
1520  u = fVertices[ipl+4]; v = fVertices[j+4];
1521  w = fVertices[(ipl+3)%4+4];
1522  }
1523  alfa = G4UniformRand();
1524  beta = G4UniformRand();
1525  lambda1=alfa*beta;
1526  lambda0=alfa-lambda1;
1527  v = v-u;
1528  w = w-u;
1529  v = u+lambda0*v+lambda1*w;
1530  }
1531  else // Lateral Plane Twisted or Not
1532  {
1533  if (chose < Surface0+Surface1) { ipl=0; }
1534  else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
1535  else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
1536  else { ipl=3; }
1537  j = (ipl+1)%4;
1538  zp = -fDz+G4UniformRand()*2*fDz;
1539  cf = 0.5*(fDz-zp)/fDz;
1540  u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
1541  v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
1542  v = u+(v-u)*G4UniformRand();
1543  }
1544  point=G4ThreeVector(v.x(),v.y(),zp);
1545 
1546  return point;
1547 }
1548 
1549 // --------------------------------------------------------------------
1550 
1552 {
1553  if(fCubicVolume != 0.) {;}
1555  return fCubicVolume;
1556 }
1557 
1558 // --------------------------------------------------------------------
1559 
1561 {
1562  if(fSurfaceArea != 0.) {;}
1563  else
1564  {
1565  std::vector<G4ThreeVector> vertices;
1566  for (G4int i=0; i<4;i++)
1567  {
1568  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1569  }
1570  for (G4int i=4; i<8;i++)
1571  {
1572  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1573  }
1574 
1575  // Surface Area of Planes(only estimation for twisted)
1576  //
1577  G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1578  vertices[2],vertices[3]);//-fDz plane
1579  G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1580  vertices[5],vertices[4]);// Lat plane
1581  G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1582  vertices[4],vertices[7]);// Lat plane
1583  G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1584  vertices[7],vertices[6]);// Lat plane
1585  G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1586  vertices[5],vertices[6]);// Lat plane
1587  G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1588  vertices[6],vertices[7]);// fDz plane
1589 
1590  // Total Surface Area
1591  //
1592  if(!fIsTwisted)
1593  {
1594  fSurfaceArea = fSurface0+fSurface1+fSurface2
1595  + fSurface3+fSurface4+fSurface5;
1596  }
1597  else
1598  {
1600  }
1601  }
1602  return fSurfaceArea;
1603 }
1604 
1605 // --------------------------------------------------------------------
1606 
1608  const G4ThreeVector& p1,
1609  const G4ThreeVector& p2,
1610  const G4ThreeVector& p3) const
1611 {
1612  // Auxiliary method for Get Surface Area of Face
1613 
1614  G4double aOne, aTwo;
1615  G4ThreeVector t, u, v, w, Area, normal;
1616 
1617  t = p2 - p1;
1618  u = p0 - p1;
1619  v = p2 - p3;
1620  w = p0 - p3;
1621 
1622  Area = w.cross(v);
1623  aOne = 0.5*Area.mag();
1624 
1625  Area = t.cross(u);
1626  aTwo = 0.5*Area.mag();
1627 
1628  return aOne + aTwo;
1629 }
1630 
1631 // --------------------------------------------------------------------
1632 
1634 {
1635  // Computes tangents of twist angles (angles between projections on XY plane
1636  // of corresponding -dz +dz edges).
1637 
1638  G4bool twisted = false;
1639  G4double dx1, dy1, dx2, dy2;
1640  G4int nv = fgkNofVertices/2;
1641 
1642  for ( G4int i=0; i<4; i++ )
1643  {
1644  dx1 = fVertices[(i+1)%nv].x()-fVertices[i].x();
1645  dy1 = fVertices[(i+1)%nv].y()-fVertices[i].y();
1646  if ( (dx1 == 0) && (dy1 == 0) ) { continue; }
1647 
1648  dx2 = fVertices[nv+(i+1)%nv].x()-fVertices[nv+i].x();
1649  dy2 = fVertices[nv+(i+1)%nv].y()-fVertices[nv+i].y();
1650 
1651  if ( dx2 == 0 && dy2 == 0 ) { continue; }
1652  G4double twist_angle = std::fabs(dy1*dx2 - dx1*dy2);
1653  if ( twist_angle < fgkTolerance ) { continue; }
1654  twisted = true;
1655  SetTwistAngle(i,twist_angle);
1656 
1657  // Check on big angles, potentially navigation problem
1658 
1659  twist_angle = std::acos( (dx1*dx2 + dy1*dy2)
1660  / (std::sqrt(dx1*dx1+dy1*dy1)
1661  * std::sqrt(dx2*dx2+dy2*dy2)) );
1662 
1663  if ( std::fabs(twist_angle) > 0.5*pi+kCarTolerance )
1664  {
1665  std::ostringstream message;
1666  message << "Twisted Angle is bigger than 90 degrees - " << GetName()
1667  << G4endl
1668  << " Potential problem of malformed Solid !" << G4endl
1669  << " TwistANGLE = " << twist_angle
1670  << "*rad for lateral plane N= " << i;
1671  G4Exception("G4GenericTrap::ComputeIsTwisted()", "GeomSolids1002",
1672  JustWarning, message);
1673  }
1674  }
1675 
1676  return twisted;
1677 }
1678 
1679 // --------------------------------------------------------------------
1680 
1681 G4bool G4GenericTrap::CheckOrder(const std::vector<G4TwoVector>& vertices) const
1682 {
1683  // Test if the vertices are in a clockwise order, if not reorder them.
1684  // Also test if they're well defined without crossing opposite segments
1685 
1686  G4bool clockwise_order=true;
1687  G4double sum1 = 0.;
1688  G4double sum2 = 0.;
1689  G4int j;
1690 
1691  for (G4int i=0; i<4; i++)
1692  {
1693  j = (i+1)%4;
1694  sum1 += vertices[i].x()*vertices[j].y() - vertices[j].x()*vertices[i].y();
1695  sum2 += vertices[i+4].x()*vertices[j+4].y()
1696  - vertices[j+4].x()*vertices[i+4].y();
1697  }
1698  if (sum1*sum2 < -fgkTolerance)
1699  {
1700  std::ostringstream message;
1701  message << "Lower/upper faces defined with opposite clockwise - "
1702  << GetName();
1703  G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids0002",
1704  FatalException, message);
1705  }
1706 
1707  if ((sum1 > 0.)||(sum2 > 0.))
1708  {
1709  std::ostringstream message;
1710  message << "Vertices must be defined in clockwise XY planes - "
1711  << GetName();
1712  G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids1001",
1713  JustWarning,message, "Re-ordering...");
1714  clockwise_order = false;
1715  }
1716 
1717  // Check for illegal crossings
1718  //
1719  G4bool illegal_cross = false;
1720  illegal_cross = IsSegCrossingZ(vertices[0],vertices[4],
1721  vertices[1],vertices[5]);
1722 
1723  if (!illegal_cross)
1724  {
1725  illegal_cross = IsSegCrossingZ(vertices[2],vertices[6],
1726  vertices[3],vertices[7]);
1727  }
1728  // +/- dZ planes
1729  if (!illegal_cross)
1730  {
1731  illegal_cross = IsSegCrossing(vertices[0],vertices[1],
1732  vertices[2],vertices[3]);
1733  }
1734  if (!illegal_cross)
1735  {
1736  illegal_cross = IsSegCrossing(vertices[0],vertices[3],
1737  vertices[1],vertices[2]);
1738  }
1739  if (!illegal_cross)
1740  {
1741  illegal_cross = IsSegCrossing(vertices[4],vertices[5],
1742  vertices[6],vertices[7]);
1743  }
1744  if (!illegal_cross)
1745  {
1746  illegal_cross = IsSegCrossing(vertices[4],vertices[7],
1747  vertices[5],vertices[6]);
1748  }
1749 
1750  if (illegal_cross)
1751  {
1752  std::ostringstream message;
1753  message << "Malformed polygone with opposite sides - " << GetName();
1754  G4Exception("G4GenericTrap::CheckOrderAndSetup()",
1755  "GeomSolids0002", FatalException, message);
1756  }
1757  return clockwise_order;
1758 }
1759 
1760 // --------------------------------------------------------------------
1761 
1762 void G4GenericTrap::ReorderVertices(std::vector<G4ThreeVector>& vertices) const
1763 {
1764  // Reorder the vector of vertices
1765 
1766  std::vector<G4ThreeVector> oldVertices(vertices);
1767 
1768  for ( G4int i=0; i < G4int(oldVertices.size()); ++i )
1769  {
1770  vertices[i] = oldVertices[oldVertices.size()-1-i];
1771  }
1772 }
1773 
1774 // --------------------------------------------------------------------
1775 
1776 G4bool
1778  const G4TwoVector& c, const G4TwoVector& d) const
1779 {
1780  // Check if segments [A,B] and [C,D] are crossing
1781 
1782  G4bool stand1 = false;
1783  G4bool stand2 = false;
1784  G4double dx1,dx2,xm=0.,ym=0.,a1=0.,a2=0.,b1=0.,b2=0.;
1785  dx1=(b-a).x();
1786  dx2=(d-c).x();
1787 
1788  if( std::fabs(dx1) < fgkTolerance ) { stand1 = true; }
1789  if( std::fabs(dx2) < fgkTolerance ) { stand2 = true; }
1790  if (!stand1)
1791  {
1792  a1 = (b.x()*a.y()-a.x()*b.y())/dx1;
1793  b1 = (b-a).y()/dx1;
1794  }
1795  if (!stand2)
1796  {
1797  a2 = (d.x()*c.y()-c.x()*d.y())/dx2;
1798  b2 = (d-c).y()/dx2;
1799  }
1800  if (stand1 && stand2)
1801  {
1802  // Segments parallel and vertical
1803  //
1804  if (std::fabs(a.x()-c.x())<fgkTolerance)
1805  {
1806  // Check if segments are overlapping
1807  //
1808  if ( ((c.y()-a.y())*(c.y()-b.y())<-fgkTolerance)
1809  || ((d.y()-a.y())*(d.y()-b.y())<-fgkTolerance)
1810  || ((a.y()-c.y())*(a.y()-d.y())<-fgkTolerance)
1811  || ((b.y()-c.y())*(b.y()-d.y())<-fgkTolerance) ) { return true; }
1812 
1813  return false;
1814  }
1815  // Different x values
1816  //
1817  return false;
1818  }
1819 
1820  if (stand1) // First segment vertical
1821  {
1822  xm = a.x();
1823  ym = a2+b2*xm;
1824  }
1825  else
1826  {
1827  if (stand2) // Second segment vertical
1828  {
1829  xm = c.x();
1830  ym = a1+b1*xm;
1831  }
1832  else // Normal crossing
1833  {
1834  if (std::fabs(b1-b2) < fgkTolerance)
1835  {
1836  // Parallel segments, are they aligned
1837  //
1838  if (std::fabs(c.y()-(a1+b1*c.x())) > fgkTolerance) { return false; }
1839 
1840  // Aligned segments, are they overlapping
1841  //
1842  if ( ((c.x()-a.x())*(c.x()-b.x())<-fgkTolerance)
1843  || ((d.x()-a.x())*(d.x()-b.x())<-fgkTolerance)
1844  || ((a.x()-c.x())*(a.x()-d.x())<-fgkTolerance)
1845  || ((b.x()-c.x())*(b.x()-d.x())<-fgkTolerance) ) { return true; }
1846 
1847  return false;
1848  }
1849  xm = (a1-a2)/(b2-b1);
1850  ym = (a1*b2-a2*b1)/(b2-b1);
1851  }
1852  }
1853 
1854  // Check if crossing point is both between A,B and C,D
1855  //
1856  G4double check = (xm-a.x())*(xm-b.x())+(ym-a.y())*(ym-b.y());
1857  if (check > -fgkTolerance) { return false; }
1858  check = (xm-c.x())*(xm-d.x())+(ym-c.y())*(ym-d.y());
1859  if (check > -fgkTolerance) { return false; }
1860 
1861  return true;
1862 }
1863 
1864 // --------------------------------------------------------------------
1865 
1866 G4bool
1868  const G4TwoVector& c, const G4TwoVector& d) const
1869 {
1870  // Check if segments [A,B] and [C,D] are crossing when
1871  // A and C are on -dZ and B and D are on +dZ
1872 
1873  // Calculate the Intersection point between two lines in 3D
1874  //
1875  G4ThreeVector temp1,temp2;
1876  G4ThreeVector v1,v2,p1,p2,p3,p4,dv;
1877  G4double q,det;
1878  p1=G4ThreeVector(a.x(),a.y(),-fDz);
1879  p2=G4ThreeVector(c.x(),c.y(),-fDz);
1880  p3=G4ThreeVector(b.x(),b.y(),fDz);
1881  p4=G4ThreeVector(d.x(),d.y(),fDz);
1882  v1=p3-p1;
1883  v2=p4-p2;
1884  dv=p2-p1;
1885 
1886  // In case of Collapsed Vertices No crossing
1887  //
1888  if( (std::fabs(dv.x()) < kCarTolerance )&&
1889  (std::fabs(dv.y()) < kCarTolerance ) ) { return false; }
1890 
1891  if( (std::fabs((p4-p3).x()) < kCarTolerance )&&
1892  (std::fabs((p4-p3).y()) < kCarTolerance ) ) { return false; }
1893 
1894  // First estimate if Intersection is possible( if det is 0)
1895  //
1896  det = dv.x()*v1.y()*v2.z()+dv.y()*v1.z()*v2.x()
1897  - dv.x()*v1.z()*v2.y()-dv.y()*v1.x()*v2.z();
1898 
1899  if (std::fabs(det)<kCarTolerance) //Intersection
1900  {
1901  temp1 = v1.cross(v2);
1902  temp2 = (p2-p1).cross(v2);
1903  if (temp1.dot(temp2) < 0) { return false; } // intersection negative
1904  q = temp1.mag();
1905 
1906  if ( q < kCarTolerance ) { return false; } // parallel lines
1907  q = ((dv).cross(v2)).mag()/q;
1908 
1909  if(q < 1.-kCarTolerance) { return true; }
1910  }
1911  return false;
1912 }
1913 
1914 // --------------------------------------------------------------------
1915 
1916 G4VFacet*
1917 G4GenericTrap::MakeDownFacet(const std::vector<G4ThreeVector>& fromVertices,
1918  G4int ind1, G4int ind2, G4int ind3) const
1919 {
1920  // Create a triangular facet from the polygon points given by indices
1921  // forming the down side ( the normal goes in -z)
1922  // Do not create facet if 2 vertices are the same
1923 
1924  if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1925  (fromVertices[ind2] == fromVertices[ind3]) ||
1926  (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
1927 
1928  std::vector<G4ThreeVector> vertices;
1929  vertices.push_back(fromVertices[ind1]);
1930  vertices.push_back(fromVertices[ind2]);
1931  vertices.push_back(fromVertices[ind3]);
1932 
1933  // first vertex most left
1934  //
1935  G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1936 
1937  if ( cross.z() > 0.0 )
1938  {
1939  // Should not happen, as vertices should have been reordered at this stage
1940 
1941  std::ostringstream message;
1942  message << "Vertices in wrong order - " << GetName();
1943  G4Exception("G4GenericTrap::MakeDownFacet", "GeomSolids0002",
1944  FatalException, message);
1945  }
1946 
1947  return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
1948 }
1949 
1950 // --------------------------------------------------------------------
1951 
1952 G4VFacet*
1953 G4GenericTrap::MakeUpFacet(const std::vector<G4ThreeVector>& fromVertices,
1954  G4int ind1, G4int ind2, G4int ind3) const
1955 {
1956  // Create a triangular facet from the polygon points given by indices
1957  // forming the upper side ( z>0 )
1958 
1959  // Do not create facet if 2 vertices are the same
1960  //
1961  if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1962  (fromVertices[ind2] == fromVertices[ind3]) ||
1963  (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
1964 
1965  std::vector<G4ThreeVector> vertices;
1966  vertices.push_back(fromVertices[ind1]);
1967  vertices.push_back(fromVertices[ind2]);
1968  vertices.push_back(fromVertices[ind3]);
1969 
1970  // First vertex most left
1971  //
1972  G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1973 
1974  if ( cross.z() < 0.0 )
1975  {
1976  // Should not happen, as vertices should have been reordered at this stage
1977 
1978  std::ostringstream message;
1979  message << "Vertices in wrong order - " << GetName();
1980  G4Exception("G4GenericTrap::MakeUpFacet", "GeomSolids0002",
1981  FatalException, message);
1982  }
1983 
1984  return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
1985 }
1986 
1987 // --------------------------------------------------------------------
1988 
1989 G4VFacet*
1991  const G4ThreeVector& downVertex1,
1992  const G4ThreeVector& upVertex1,
1993  const G4ThreeVector& upVertex0) const
1994 {
1995  // Creates a triangular facet from the polygon points given by indices
1996  // forming the upper side ( z>0 )
1997 
1998  if ( (downVertex0 == downVertex1) && (upVertex0 == upVertex1) )
1999  {
2000  return 0;
2001  }
2002 
2003  if ( downVertex0 == downVertex1 )
2004  {
2005  return new G4TriangularFacet(downVertex0, upVertex1, upVertex0, ABSOLUTE);
2006  }
2007 
2008  if ( upVertex0 == upVertex1 )
2009  {
2010  return new G4TriangularFacet(downVertex0, downVertex1, upVertex0, ABSOLUTE);
2011  }
2012 
2013  return new G4QuadrangularFacet(downVertex0, downVertex1,
2014  upVertex1, upVertex0, ABSOLUTE);
2015 }
2016 
2017 // --------------------------------------------------------------------
2018 
2020 {
2021  // 3D vertices
2022  //
2023  G4int nv = fgkNofVertices/2;
2024  std::vector<G4ThreeVector> downVertices;
2025  for ( G4int i=0; i<nv; i++ )
2026  {
2027  downVertices.push_back(G4ThreeVector(fVertices[i].x(),
2028  fVertices[i].y(), -fDz));
2029  }
2030 
2031  std::vector<G4ThreeVector> upVertices;
2032  for ( G4int i=nv; i<2*nv; i++ )
2033  {
2034  upVertices.push_back(G4ThreeVector(fVertices[i].x(),
2035  fVertices[i].y(), fDz));
2036  }
2037 
2038  // Reorder vertices if they are not ordered anti-clock wise
2039  //
2040  G4ThreeVector cross
2041  = (downVertices[1]-downVertices[0]).cross(downVertices[2]-downVertices[1]);
2042  G4ThreeVector cross1
2043  = (upVertices[1]-upVertices[0]).cross(upVertices[2]-upVertices[1]);
2044  if ( (cross.z() > 0.0) || (cross1.z() > 0.0) )
2045  {
2046  ReorderVertices(downVertices);
2047  ReorderVertices(upVertices);
2048  }
2049 
2050  G4TessellatedSolid* tessellatedSolid = new G4TessellatedSolid(GetName());
2051 
2052  G4VFacet* facet = 0;
2053  facet = MakeDownFacet(downVertices, 0, 1, 2);
2054  if (facet) { tessellatedSolid->AddFacet( facet ); }
2055  facet = MakeDownFacet(downVertices, 0, 2, 3);
2056  if (facet) { tessellatedSolid->AddFacet( facet ); }
2057  facet = MakeUpFacet(upVertices, 0, 2, 1);
2058  if (facet) { tessellatedSolid->AddFacet( facet ); }
2059  facet = MakeUpFacet(upVertices, 0, 3, 2);
2060  if (facet) { tessellatedSolid->AddFacet( facet ); }
2061 
2062  // The quadrangular sides
2063  //
2064  for ( G4int i = 0; i < nv; ++i )
2065  {
2066  G4int j = (i+1) % nv;
2067  facet = MakeSideFacet(downVertices[j], downVertices[i],
2068  upVertices[i], upVertices[j]);
2069 
2070  if ( facet ) { tessellatedSolid->AddFacet( facet ); }
2071  }
2072 
2073  tessellatedSolid->SetSolidClosed(true);
2074 
2075  return tessellatedSolid;
2076 }
2077 
2078 // --------------------------------------------------------------------
2079 
2081 {
2082  // Computes bounding box for a shape.
2083 
2084  G4double minX, maxX, minY, maxY;
2085  minX = maxX = fVertices[0].x();
2086  minY = maxY = fVertices[0].y();
2087 
2088  for (G4int i=1; i< fgkNofVertices; i++)
2089  {
2090  if (minX>fVertices[i].x()) { minX=fVertices[i].x(); }
2091  if (maxX<fVertices[i].x()) { maxX=fVertices[i].x(); }
2092  if (minY>fVertices[i].y()) { minY=fVertices[i].y(); }
2093  if (maxY<fVertices[i].y()) { maxY=fVertices[i].y(); }
2094  }
2095  fMinBBoxVector = G4ThreeVector(minX,minY,-fDz);
2096  fMaxBBoxVector = G4ThreeVector(maxX,maxY, fDz);
2097 }
2098 
2099 // --------------------------------------------------------------------
2100 
2102 {
2103 
2104 #ifdef G4TESS_TEST
2105  if ( fTessellatedSolid )
2106  {
2107  return fTessellatedSolid->GetPolyhedron();
2108  }
2109 #endif
2110 
2111  if ( (!fpPolyhedron)
2114  fpPolyhedron->GetNumberOfRotationSteps()) )
2115  {
2116  G4AutoLock l(&polyhedronMutex);
2117  delete fpPolyhedron;
2119  fRebuildPolyhedron = false;
2120  l.unlock();
2121  }
2122  return fpPolyhedron;
2123 }
2124 
2125 // --------------------------------------------------------------------
2126 
2128 {
2129 
2130 #ifdef G4TESS_TEST
2131  if ( fTessellatedSolid )
2132  {
2133  return fTessellatedSolid->DescribeYourselfTo(scene);
2134  }
2135 #endif
2136 
2137  scene.AddSolid(*this);
2138 }
2139 
2140 // --------------------------------------------------------------------
2141 
2143 {
2144  // Computes bounding vectors for the shape
2145 
2146 #ifdef G4TESS_TEST
2147  if ( fTessellatedSolid )
2148  {
2149  return fTessellatedSolid->GetExtent();
2150  }
2151 #endif
2152 
2153  G4double Dx,Dy;
2154  G4ThreeVector minVec = GetMinimumBBox();
2155  G4ThreeVector maxVec = GetMaximumBBox();
2156  Dx = 0.5*(maxVec.x()- minVec.x());
2157  Dy = 0.5*(maxVec.y()- minVec.y());
2158 
2159  return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz);
2160 }
2161 
2162 // --------------------------------------------------------------------
2163 
2165 {
2166 
2167 #ifdef G4TESS_TEST
2168  if ( fTessellatedSolid )
2169  {
2171  }
2172 #endif
2173 
2174  // Approximation of Twisted Side
2175  // Construct extra Points, if Twisted Side
2176  //
2177  G4PolyhedronArbitrary* polyhedron;
2178  size_t nVertices, nFacets;
2179 
2180  G4int subdivisions=0;
2181  G4int i;
2182  if(fIsTwisted)
2183  {
2184  if ( GetVisSubdivisions()!= 0 )
2185  {
2186  subdivisions=GetVisSubdivisions();
2187  }
2188  else
2189  {
2190  // Estimation of Number of Subdivisions for smooth visualisation
2191  //
2192  G4double maxTwist=0.;
2193  for(i=0; i<4; i++)
2194  {
2195  if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
2196  }
2197 
2198  // Computes bounding vectors for the shape
2199  //
2200  G4double Dx,Dy;
2201  G4ThreeVector minVec = GetMinimumBBox();
2202  G4ThreeVector maxVec = GetMaximumBBox();
2203  Dx = 0.5*(maxVec.x()- minVec.y());
2204  Dy = 0.5*(maxVec.y()- minVec.y());
2205  if (Dy > Dx) { Dx=Dy; }
2206 
2207  subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
2208  if (subdivisions<4) { subdivisions=4; }
2209  if (subdivisions>30) { subdivisions=30; }
2210  }
2211  }
2212  G4int sub4=4*subdivisions;
2213  nVertices = 8+subdivisions*4;
2214  nFacets = 6+subdivisions*4;
2215  G4double cf=1./(subdivisions+1);
2216  polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
2217 
2218  // Add Vertex
2219  //
2220  for (i=0;i<4;i++)
2221  {
2222  polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
2223  fVertices[i].y(),-fDz));
2224  }
2225  for( i=0;i<subdivisions;i++)
2226  {
2227  for(G4int j=0;j<4;j++)
2228  {
2229  G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
2230  polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
2231  }
2232  }
2233  for (i=4;i<8;i++)
2234  {
2235  polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
2236  fVertices[i].y(),fDz));
2237  }
2238 
2239  // Add Facets
2240  //
2241  polyhedron->AddFacet(1,4,3,2); //Z-plane
2242  for (i=0;i<subdivisions+1;i++)
2243  {
2244  G4int is=i*4;
2245  polyhedron->AddFacet(5+is,8+is,4+is,1+is);
2246  polyhedron->AddFacet(8+is,7+is,3+is,4+is);
2247  polyhedron->AddFacet(7+is,6+is,2+is,3+is);
2248  polyhedron->AddFacet(6+is,5+is,1+is,2+is);
2249  }
2250  polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4); //Z-plane
2251 
2252  polyhedron->SetReferences();
2253  polyhedron->InvertFacets();
2254 
2255  return (G4Polyhedron*) polyhedron;
2256 }
2257 
2258 // --------------------------------------------------------------------
2259 
2260 #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