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