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