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