source: trunk/source/geometry/solids/specific/src/G4GenericTrap.cc@ 1347

Last change on this file since 1347 was 1337, checked in by garnier, 15 years ago

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