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