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