Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4Trap.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// Implementation for G4Trap class
27//
28// 21.03.95 P.Kent: Modified for `tolerant' geometry
29// 09.09.96 V.Grichine: Final modifications before to commit
30// 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters
31// 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal
32// 18.04.17 E.Tcherniaev: complete revision, speed-up
33// --------------------------------------------------------------------
34
35#include "G4Trap.hh"
36
37#if !defined(G4GEOM_USE_UTRAP)
38
39#include "globals.hh"
40#include "G4GeomTools.hh"
41
42#include "G4VoxelLimits.hh"
43#include "G4AffineTransform.hh"
44#include "G4BoundingEnvelope.hh"
45
46#include "G4VPVParameterisation.hh"
47
48#include "G4QuickRand.hh"
49
50#include "G4VGraphicsScene.hh"
51#include "G4Polyhedron.hh"
52#include "G4AutoLock.hh"
53
54namespace
55{
56 G4Mutex trapMutex = G4MUTEX_INITIALIZER;
57}
58
59using namespace CLHEP;
60
61//////////////////////////////////////////////////////////////////////////
62//
63// Constructor - check and set half-widths as well as angles:
64// final check of coplanarity
65
66G4Trap::G4Trap( const G4String& pName,
67 G4double pDz,
68 G4double pTheta, G4double pPhi,
69 G4double pDy1, G4double pDx1, G4double pDx2,
70 G4double pAlp1,
71 G4double pDy2, G4double pDx3, G4double pDx4,
72 G4double pAlp2 )
73 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
74{
75 fDz = pDz;
76 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
77 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
78
79 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
80 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
81
82 CheckParameters();
83 MakePlanes();
84}
85
86//////////////////////////////////////////////////////////////////////////
87//
88// Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
89// which are its vertices. Checking of planarity with preparation of
90// fPlanes[] and than calculation of other members
91
92G4Trap::G4Trap( const G4String& pName,
93 const G4ThreeVector pt[8] )
94 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
95{
96 // Start with check of centering - the center of gravity trap line
97 // should cross the origin of frame
98 //
99 if ( pt[0].z() >= 0
100 || pt[0].z() != pt[1].z()
101 || pt[0].z() != pt[2].z()
102 || pt[0].z() != pt[3].z()
103
104 || pt[4].z() <= 0
105 || pt[4].z() != pt[5].z()
106 || pt[4].z() != pt[6].z()
107 || pt[4].z() != pt[7].z()
108
109 || std::fabs( pt[0].z() + pt[4].z() ) >= kCarTolerance
110
111 || pt[0].y() != pt[1].y()
112 || pt[2].y() != pt[3].y()
113 || pt[4].y() != pt[5].y()
114 || pt[6].y() != pt[7].y()
115
116 || std::fabs(pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y()) >= kCarTolerance
117 || std::fabs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() +
118 pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) >= kCarTolerance )
119 {
120 std::ostringstream message;
121 message << "Invalid vertice coordinates for Solid: " << GetName();
122 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
123 FatalException, message);
124 }
125
126 // Set parameters
127 //
128 fDz = (pt[7]).z();
129
130 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
131 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
132 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
133 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
134
135 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
136 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
137 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
138 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
139
140 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
141 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
142
143 CheckParameters();
144 MakePlanes(pt);
145}
146
147//////////////////////////////////////////////////////////////////////////
148//
149// Constructor for Right Angular Wedge from STEP
150
151G4Trap::G4Trap( const G4String& pName,
152 G4double pZ,
153 G4double pY,
154 G4double pX, G4double pLTX )
155 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
156{
157 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi = 0;
158 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLTX; fTalpha1 = 0.5*(pLTX - pX)/pY;
159 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; fTalpha2 = fTalpha1;
160
161 CheckParameters();
162 MakePlanes();
163}
164
165//////////////////////////////////////////////////////////////////////////
166//
167// Constructor for G4Trd
168
169G4Trap::G4Trap( const G4String& pName,
170 G4double pDx1, G4double pDx2,
171 G4double pDy1, G4double pDy2,
172 G4double pDz )
173 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance), fTrapType(0)
174{
175 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = 0;
176 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalpha1 = 0;
177 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalpha2 = 0;
178
179 CheckParameters();
180 MakePlanes();
181}
182
183//////////////////////////////////////////////////////////////////////////
184//
185// Constructor for G4Para
186
187G4Trap::G4Trap( const G4String& pName,
188 G4double pDx, G4double pDy,
189 G4double pDz,
190 G4double pAlpha,
191 G4double pTheta, G4double pPhi )
192 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
193{
194 fDz = pDz;
195 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
196 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
197
198 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 = std::tan(pAlpha);
199 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 = fTalpha1;
200
201 CheckParameters();
202 MakePlanes();
203}
204
205//////////////////////////////////////////////////////////////////////////
206//
207// Nominal constructor for G4Trap whose parameters are to be set by
208// a G4VParamaterisation later. Check and set half-widths as well as
209// angles: final check of coplanarity
210
211G4Trap::G4Trap( const G4String& pName )
212 : G4CSGSolid (pName), halfCarTolerance(0.5*kCarTolerance),
213 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
214 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
215 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
216{
217 MakePlanes();
218}
219
220//////////////////////////////////////////////////////////////////////////
221//
222// Fake default constructor - sets only member data and allocates memory
223// for usage restricted to object persistency.
224//
225G4Trap::G4Trap( __void__& a )
226 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance),
227 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
228 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
229 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
230{
231 MakePlanes();
232}
233
234//////////////////////////////////////////////////////////////////////////
235//
236// Copy constructor
237
238G4Trap::G4Trap(const G4Trap& rhs)
239 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
243{
244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; }
246 fTrapType = rhs.fTrapType;
247}
248
249//////////////////////////////////////////////////////////////////////////
250//
251// Assignment operator
252
253G4Trap& G4Trap::operator = (const G4Trap& rhs)
254{
255 // Check assignment to self
256 //
257 if (this == &rhs) { return *this; }
258
259 // Copy base class data
260 //
261 G4CSGSolid::operator=(rhs);
262
263 // Copy data
264 //
265 halfCarTolerance = rhs.halfCarTolerance;
266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; }
271 fTrapType = rhs.fTrapType;
272 return *this;
273}
274
275//////////////////////////////////////////////////////////////////////////
276//
277// Set all parameters, as for constructor - check and set half-widths
278// as well as angles: final check of coplanarity
279
280void G4Trap::SetAllParameters ( G4double pDz,
281 G4double pTheta,
282 G4double pPhi,
283 G4double pDy1,
284 G4double pDx1,
285 G4double pDx2,
286 G4double pAlp1,
287 G4double pDy2,
288 G4double pDx3,
289 G4double pDx4,
290 G4double pAlp2 )
291{
292 // Reset data of the base class
293 fCubicVolume = 0;
294 fSurfaceArea = 0;
295 fRebuildPolyhedron = true;
296
297 // Set parameters
298 fDz = pDz;
299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
301
302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
304
305 CheckParameters();
306 MakePlanes();
307}
308
309//////////////////////////////////////////////////////////////////////////
310//
311// Check length parameters
312
313void G4Trap::CheckParameters()
314{
315 if (fDz<=0 ||
316 fDy1<=0 || fDx1<=0 || fDx2<=0 ||
317 fDy2<=0 || fDx3<=0 || fDx4<=0)
318 {
319 std::ostringstream message;
320 message << "Invalid Length Parameters for Solid: " << GetName()
321 << "\n X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4
322 << "\n Y - " <<fDy1<<", "<<fDy2
323 << "\n Z - " <<fDz;
324 G4Exception("G4Trap::CheckParameters()", "GeomSolids0002",
325 FatalException, message);
326 }
327}
328
329//////////////////////////////////////////////////////////////////////////
330//
331// Compute vertices and set side planes
332
333void G4Trap::MakePlanes()
334{
335 G4double DzTthetaCphi = fDz*fTthetaCphi;
336 G4double DzTthetaSphi = fDz*fTthetaSphi;
337 G4double Dy1Talpha1 = fDy1*fTalpha1;
338 G4double Dy2Talpha2 = fDy2*fTalpha2;
339
340 G4ThreeVector pt[8] =
341 {
342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz),
343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz),
344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz),
345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz),
346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz),
347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz),
348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz),
349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz)
350 };
351
352 MakePlanes(pt);
353}
354
355//////////////////////////////////////////////////////////////////////////
356//
357// Set side planes, check planarity
358
359void G4Trap::MakePlanes(const G4ThreeVector pt[8])
360{
361 constexpr G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} };
362 const static G4String side[4] = { "~-Y", "~+Y", "~-X", "~+X" };
363
364 for (G4int i=0; i<4; ++i)
365 {
366 if (MakePlane(pt[iface[i][0]],
367 pt[iface[i][1]],
368 pt[iface[i][2]],
369 pt[iface[i][3]],
370 fPlanes[i])) { continue; }
371
372 // Non planar side face
373 G4ThreeVector normal(fPlanes[i].a,fPlanes[i].b,fPlanes[i].c);
374 G4double dmax = 0;
375 for (G4int k=0; k<4; ++k)
376 {
377 G4double dist = normal.dot(pt[iface[i][k]]) + fPlanes[i].d;
378 if (std::abs(dist) > std::abs(dmax)) { dmax = dist; }
379 }
380 std::ostringstream message;
381 message << "Side face " << side[i] << " is not planar for solid: "
382 << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n";
383 StreamInfo(message);
384 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
385 FatalException, message);
386 }
387
388 // Re-compute parameters
389 SetCachedValues();
390}
391
392//////////////////////////////////////////////////////////////////////////
393//
394// Calculate the coef's of the plane p1->p2->p3->p4->p1
395// where the ThreeVectors 1-4 are in anti-clockwise order when viewed
396// from infront of the plane (i.e. from normal direction).
397//
398// Return true if the points are coplanar, false otherwise
399
400G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401 const G4ThreeVector& p2,
402 const G4ThreeVector& p3,
403 const G4ThreeVector& p4,
404 TrapSidePlane& plane )
405{
406 G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit();
407 if (std::abs(normal.x()) < DBL_EPSILON) { normal.setX(0); }
408 if (std::abs(normal.y()) < DBL_EPSILON) { normal.setY(0); }
409 if (std::abs(normal.z()) < DBL_EPSILON) { normal.setZ(0); }
410 normal = normal.unit();
411
412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25;
413 plane.a = normal.x();
414 plane.b = normal.y();
415 plane.c = normal.z();
416 plane.d = -normal.dot(centre);
417
418 // compute distances and check planarity
419 G4double d1 = std::abs(normal.dot(p1) + plane.d);
420 G4double d2 = std::abs(normal.dot(p2) + plane.d);
421 G4double d3 = std::abs(normal.dot(p3) + plane.d);
422 G4double d4 = std::abs(normal.dot(p4) + plane.d);
423 G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4);
424
425 return dmax <= 1000 * kCarTolerance;
426}
427
428//////////////////////////////////////////////////////////////////////////
429//
430// Recompute parameters using planes
431
432void G4Trap::SetCachedValues()
433{
434 // Set indeces
435 constexpr G4int iface[6][4] =
436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
437
438 // Get vertices
439 G4ThreeVector pt[8];
440 GetVertices(pt);
441
442 // Set face areas
443 for (G4int i=0; i<6; ++i)
444 {
445 fAreas[i] = G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
446 pt[iface[i][1]],
447 pt[iface[i][2]],
448 pt[iface[i][3]]).mag();
449 }
450 for (G4int i=1; i<6; ++i) { fAreas[i] += fAreas[i - 1]; }
451
452 // Define type of trapezoid
453 fTrapType = 0;
454 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 &&
455 std::abs(fPlanes[0].a) < DBL_EPSILON &&
456 std::abs(fPlanes[0].c) < DBL_EPSILON &&
457 std::abs(fPlanes[1].a) < DBL_EPSILON &&
458 std::abs(fPlanes[1].c) < DBL_EPSILON)
459 {
460 fTrapType = 1; // YZ section is a rectangle ...
461 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
462 std::abs(fPlanes[2].c - fPlanes[3].c) < DBL_EPSILON &&
463 fPlanes[2].b == 0 &&
464 fPlanes[3].b == 0)
465 {
466 fTrapType = 2; // ... and XZ section is a isosceles trapezoid
467 fPlanes[2].a = -fPlanes[3].a;
468 fPlanes[2].c = fPlanes[3].c;
469 }
470 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
471 std::abs(fPlanes[2].b - fPlanes[3].b) < DBL_EPSILON &&
472 fPlanes[2].c == 0 &&
473 fPlanes[3].c == 0)
474 {
475 fTrapType = 3; // ... and XY section is a isosceles trapezoid
476 fPlanes[2].a = -fPlanes[3].a;
477 fPlanes[2].b = fPlanes[3].b;
478 }
479 }
480}
481
482//////////////////////////////////////////////////////////////////////////
483//
484// Dispatch to parameterisation for replication mechanism dimension
485// computation & modification.
486
487void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
488 const G4int n,
489 const G4VPhysicalVolume* pRep )
490{
491 p->ComputeDimensions(*this,n,pRep);
492}
493
494//////////////////////////////////////////////////////////////////////////
495//
496// Get bounding box
497
498void G4Trap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
499{
500 G4ThreeVector pt[8];
501 GetVertices(pt);
502
503 G4double xmin = kInfinity, xmax = -kInfinity;
504 G4double ymin = kInfinity, ymax = -kInfinity;
505 for (const auto & i : pt)
506 {
507 G4double x = i.x();
508 if (x < xmin) { xmin = x; }
509 if (x > xmax) { xmax = x; }
510 G4double y = i.y();
511 if (y < ymin) { ymin = y; }
512 if (y > ymax) { ymax = y; }
513 }
514
515 G4double dz = GetZHalfLength();
516 pMin.set(xmin,ymin,-dz);
517 pMax.set(xmax,ymax, dz);
518
519 // Check correctness of the bounding box
520 //
521 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
522 {
523 std::ostringstream message;
524 message << "Bad bounding box (min >= max) for solid: "
525 << GetName() << " !"
526 << "\npMin = " << pMin
527 << "\npMax = " << pMax;
528 G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001",
529 JustWarning, message);
530 DumpInfo();
531 }
532}
533
534//////////////////////////////////////////////////////////////////////////
535//
536// Calculate extent under transform and specified limit
537
538G4bool G4Trap::CalculateExtent( const EAxis pAxis,
539 const G4VoxelLimits& pVoxelLimit,
540 const G4AffineTransform& pTransform,
541 G4double& pMin, G4double& pMax) const
542{
543 G4ThreeVector bmin, bmax;
544 G4bool exist;
545
546 // Check bounding box (bbox)
547 //
548 BoundingLimits(bmin,bmax);
549 G4BoundingEnvelope bbox(bmin,bmax);
550#ifdef G4BBOX_EXTENT
551 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
552#endif
553 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
554 {
555 return exist = pMin < pMax;
556 }
557
558 // Set bounding envelope (benv) and calculate extent
559 //
560 G4ThreeVector pt[8];
561 GetVertices(pt);
562
563 G4ThreeVectorList baseA(4), baseB(4);
564 baseA[0] = pt[0];
565 baseA[1] = pt[1];
566 baseA[2] = pt[3];
567 baseA[3] = pt[2];
568
569 baseB[0] = pt[4];
570 baseB[1] = pt[5];
571 baseB[2] = pt[7];
572 baseB[3] = pt[6];
573
574 std::vector<const G4ThreeVectorList *> polygons(2);
575 polygons[0] = &baseA;
576 polygons[1] = &baseB;
577
578 G4BoundingEnvelope benv(bmin,bmax,polygons);
579 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
580 return exist;
581}
582
583//////////////////////////////////////////////////////////////////////////
584//
585// Return whether point is inside/outside/on_surface
586
587EInside G4Trap::Inside( const G4ThreeVector& p ) const
588{
589 switch (fTrapType)
590 {
591 case 0: // General case
592 {
593 G4double dz = std::abs(p.z())-fDz;
594 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
595 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
596 G4double dy = std::max(dz,std::max(dy1,dy2));
597
598 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
599 + fPlanes[2].c*p.z()+fPlanes[2].d;
600 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
601 + fPlanes[3].c*p.z()+fPlanes[3].d;
602 G4double dist = std::max(dy,std::max(dx1,dx2));
603
604 return (dist > halfCarTolerance) ? kOutside :
605 ((dist > -halfCarTolerance) ? kSurface : kInside);
606 }
607 case 1: // YZ section is a rectangle
608 {
609 G4double dz = std::abs(p.z())-fDz;
610 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
611 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
612 + fPlanes[2].c*p.z()+fPlanes[2].d;
613 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
614 + fPlanes[3].c*p.z()+fPlanes[3].d;
615 G4double dist = std::max(dy,std::max(dx1,dx2));
616
617 return (dist > halfCarTolerance) ? kOutside :
618 ((dist > -halfCarTolerance) ? kSurface : kInside);
619 }
620 case 2: // YZ section is a rectangle and
621 { // XZ section is an isosceles trapezoid
622 G4double dz = std::abs(p.z())-fDz;
623 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
624 G4double dx = fPlanes[3].a*std::abs(p.x())
625 + fPlanes[3].c*p.z()+fPlanes[3].d;
626 G4double dist = std::max(dy,dx);
627
628 return (dist > halfCarTolerance) ? kOutside :
629 ((dist > -halfCarTolerance) ? kSurface : kInside);
630 }
631 case 3: // YZ section is a rectangle and
632 { // XY section is an isosceles trapezoid
633 G4double dz = std::abs(p.z())-fDz;
634 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
635 G4double dx = fPlanes[3].a*std::abs(p.x())
636 + fPlanes[3].b*p.y()+fPlanes[3].d;
637 G4double dist = std::max(dy,dx);
638
639 return (dist > halfCarTolerance) ? kOutside :
640 ((dist > -halfCarTolerance) ? kSurface : kInside);
641 }
642 }
643 return kOutside;
644}
645
646//////////////////////////////////////////////////////////////////////////
647//
648// Determine side, and return corresponding normal
649
650G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
651{
652 G4double nx = 0, ny = 0, nz = 0;
653 G4double dz = std::abs(p.z()) - fDz;
654 nz = std::copysign(G4double(std::abs(dz) <= halfCarTolerance), p.z());
655
656 switch (fTrapType)
657 {
658 case 0: // General case
659 {
660 for (G4int i=0; i<2; ++i)
661 {
662 G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
663 if (std::abs(dy) > halfCarTolerance) { continue; }
664 ny = fPlanes[i].b;
665 nz += fPlanes[i].c;
666 break;
667 }
668 for (G4int i=2; i<4; ++i)
669 {
670 G4double dx = fPlanes[i].a*p.x() +
671 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
672 if (std::abs(dx) > halfCarTolerance) { continue; }
673 nx = fPlanes[i].a;
674 ny += fPlanes[i].b;
675 nz += fPlanes[i].c;
676 break;
677 }
678 break;
679 }
680 case 1: // YZ section - rectangle
681 {
682 G4double dy = std::abs(p.y()) + fPlanes[1].d;
683 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
684 for (G4int i=2; i<4; ++i)
685 {
686 G4double dx = fPlanes[i].a*p.x() +
687 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
688 if (std::abs(dx) > halfCarTolerance) { continue; }
689 nx = fPlanes[i].a;
690 ny += fPlanes[i].b;
691 nz += fPlanes[i].c;
692 break;
693 }
694 break;
695 }
696 case 2: // YZ section - rectangle, XZ section - isosceles trapezoid
697 {
698 G4double dy = std::abs(p.y()) + fPlanes[1].d;
699 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
700 G4double dx = fPlanes[3].a*std::abs(p.x()) +
701 fPlanes[3].c*p.z() + fPlanes[3].d;
702 G4double k = static_cast<G4double>(std::abs(dx) <= halfCarTolerance);
703 nx = std::copysign(k, p.x())*fPlanes[3].a;
704 nz += k*fPlanes[3].c;
705 break;
706 }
707 case 3: // YZ section - rectangle, XY section - isosceles trapezoid
708 {
709 G4double dy = std::abs(p.y()) + fPlanes[1].d;
710 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
711 G4double dx = fPlanes[3].a*std::abs(p.x()) +
712 fPlanes[3].b*p.y() + fPlanes[3].d;
713 G4double k = static_cast<G4double>(std::abs(dx) <= halfCarTolerance);
714 nx = std::copysign(k, p.x())*fPlanes[3].a;
715 ny += k*fPlanes[3].b;
716 break;
717 }
718 }
719
720 // Return normal
721 //
722 G4double mag2 = nx*nx + ny*ny + nz*nz;
723 if (mag2 == 1)
724 {
725 return { nx,ny,nz };
726 }
727 if (mag2 != 0)
728 {
729 return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
730 }
731
732 // Point is not on the surface
733 //
734#ifdef G4CSGDEBUG
735 std::ostringstream message;
736 G4long oldprc = message.precision(16);
737 message << "Point p is not on surface (!?) of solid: "
738 << GetName() << G4endl;
739 message << "Position:\n";
740 message << " p.x() = " << p.x()/mm << " mm\n";
741 message << " p.y() = " << p.y()/mm << " mm\n";
742 message << " p.z() = " << p.z()/mm << " mm";
743 G4cout.precision(oldprc) ;
744 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
745 JustWarning, message );
746 DumpInfo();
747#endif
748
749 return ApproxSurfaceNormal(p);
750}
751
752//////////////////////////////////////////////////////////////////////////
753//
754// Algorithm for SurfaceNormal() following the original specification
755// for points not on the surface
756
757G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
758{
759 G4double dist = -DBL_MAX;
760 G4int iside = 0;
761 for (G4int i=0; i<4; ++i)
762 {
763 G4double d = fPlanes[i].a*p.x() +
764 fPlanes[i].b*p.y() +
765 fPlanes[i].c*p.z() + fPlanes[i].d;
766 if (d > dist) { dist = d; iside = i; }
767 }
768
769 G4double distz = std::abs(p.z()) - fDz;
770 if (dist > distz)
771 {
772 return { fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c };
773 }
774
775 return { 0, 0, (G4double)((p.z() < 0) ? -1 : 1) };
776}
777
778//////////////////////////////////////////////////////////////////////////
779//
780// Calculate distance to shape from outside
781// - return kInfinity if no intersection
782
783G4double G4Trap::DistanceToIn(const G4ThreeVector& p,
784 const G4ThreeVector& v ) const
785{
786 // Z intersections
787 //
788 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
789 {
790 return kInfinity;
791 }
792 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
793 G4double dz = (invz < 0) ? fDz : -fDz;
794 G4double tzmin = (p.z() + dz)*invz;
795 G4double tzmax = (p.z() - dz)*invz;
796
797 // Y intersections
798 //
799 G4double tymin = 0, tymax = DBL_MAX;
800 G4int i = 0;
801 for ( ; i<2; ++i)
802 {
803 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
804 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
805 if (dist >= -halfCarTolerance)
806 {
807 if (cosa >= 0) { return kInfinity; }
808 G4double tmp = -dist/cosa;
809 if (tymin < tmp) { tymin = tmp; }
810 }
811 else if (cosa > 0)
812 {
813 G4double tmp = -dist/cosa;
814 if (tymax > tmp) { tymax = tmp; }
815 }
816 }
817
818 // Z intersections
819 //
820 G4double txmin = 0, txmax = DBL_MAX;
821 for ( ; i<4; ++i)
822 {
823 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
824 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() +
825 fPlanes[i].d;
826 if (dist >= -halfCarTolerance)
827 {
828 if (cosa >= 0) { return kInfinity; }
829 G4double tmp = -dist/cosa;
830 if (txmin < tmp) { txmin = tmp; }
831 }
832 else if (cosa > 0)
833 {
834 G4double tmp = -dist/cosa;
835 if (txmax > tmp) { txmax = tmp; }
836 }
837 }
838
839 // Find distance
840 //
841 G4double tmin = std::max(std::max(txmin,tymin),tzmin);
842 G4double tmax = std::min(std::min(txmax,tymax),tzmax);
843
844 if (tmax <= tmin + halfCarTolerance) { return kInfinity; } // touch or no hit
845
846 return (tmin < halfCarTolerance ) ? 0. : tmin;
847}
848
849//////////////////////////////////////////////////////////////////////////
850//
851// Calculate exact shortest distance to any boundary from outside
852// This is the best fast estimation of the shortest distance to trap
853// - return 0 if point is inside
854
855G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
856{
857 switch (fTrapType)
858 {
859 case 0: // General case
860 {
861 G4double dz = std::abs(p.z())-fDz;
862 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
863 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
864 G4double dy = std::max(dz,std::max(dy1,dy2));
865
866 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
867 + fPlanes[2].c*p.z()+fPlanes[2].d;
868 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
869 + fPlanes[3].c*p.z()+fPlanes[3].d;
870 G4double dist = std::max(dy,std::max(dx1,dx2));
871 return (dist > 0) ? dist : 0.;
872 }
873 case 1: // YZ section is a rectangle
874 {
875 G4double dz = std::abs(p.z())-fDz;
876 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
877 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
878 + fPlanes[2].c*p.z()+fPlanes[2].d;
879 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
880 + fPlanes[3].c*p.z()+fPlanes[3].d;
881 G4double dist = std::max(dy,std::max(dx1,dx2));
882 return (dist > 0) ? dist : 0.;
883 }
884 case 2: // YZ section is a rectangle and
885 { // XZ section is an isosceles trapezoid
886 G4double dz = std::abs(p.z())-fDz;
887 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
888 G4double dx = fPlanes[3].a*std::abs(p.x())
889 + fPlanes[3].c*p.z()+fPlanes[3].d;
890 G4double dist = std::max(dy,dx);
891 return (dist > 0) ? dist : 0.;
892 }
893 case 3: // YZ section is a rectangle and
894 { // XY section is an isosceles trapezoid
895 G4double dz = std::abs(p.z())-fDz;
896 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
897 G4double dx = fPlanes[3].a*std::abs(p.x())
898 + fPlanes[3].b*p.y()+fPlanes[3].d;
899 G4double dist = std::max(dy,dx);
900 return (dist > 0) ? dist : 0.;
901 }
902 }
903 return 0.;
904}
905
906//////////////////////////////////////////////////////////////////////////
907//
908// Calculate distance to surface of shape from inside and
909// find normal at exit point, if required
910// - when leaving the surface, return 0
911
912G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
913 const G4bool calcNorm,
914 G4bool* validNorm, G4ThreeVector* n) const
915{
916 // Z intersections
917 //
918 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
919 {
920 if (calcNorm)
921 {
922 *validNorm = true;
923 n->set(0, 0, (p.z() < 0) ? -1 : 1);
924 }
925 return 0;
926 }
927 G4double vz = v.z();
928 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
929 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
930
931 // Y intersections
932 //
933 G4int i = 0;
934 for ( ; i<2; ++i)
935 {
936 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
937 if (cosa > 0)
938 {
939 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
940 if (dist >= -halfCarTolerance)
941 {
942 if (calcNorm)
943 {
944 *validNorm = true;
945 n->set(0, fPlanes[i].b, fPlanes[i].c);
946 }
947 return 0;
948 }
949 G4double tmp = -dist/cosa;
950 if (tmax > tmp) { tmax = tmp; iside = i; }
951 }
952 }
953
954 // X intersections
955 //
956 for ( ; i<4; ++i)
957 {
958 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
959 if (cosa > 0)
960 {
961 G4double dist = fPlanes[i].a*p.x() +
962 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
963 if (dist >= -halfCarTolerance)
964 {
965 if (calcNorm)
966 {
967 *validNorm = true;
968 n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
969 }
970 return 0;
971 }
972 G4double tmp = -dist/cosa;
973 if (tmax > tmp) { tmax = tmp; iside = i; }
974 }
975 }
976
977 // Set normal, if required, and return distance
978 //
979 if (calcNorm)
980 {
981 *validNorm = true;
982 if (iside < 0)
983 {
984 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
985 }
986 else
987 {
988 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
989 }
990 }
991 return tmax;
992}
993
994//////////////////////////////////////////////////////////////////////////
995//
996// Calculate exact shortest distance to any boundary from inside
997// - Returns 0 is ThreeVector outside
998
999G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1000{
1001#ifdef G4CSGDEBUG
1002 if( Inside(p) == kOutside )
1003 {
1004 std::ostringstream message;
1005 G4long oldprc = message.precision(16);
1006 message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
1007 message << "Position:\n";
1008 message << " p.x() = " << p.x()/mm << " mm\n";
1009 message << " p.y() = " << p.y()/mm << " mm\n";
1010 message << " p.z() = " << p.z()/mm << " mm";
1011 G4cout.precision(oldprc);
1012 G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
1013 JustWarning, message );
1014 DumpInfo();
1015 }
1016#endif
1017 switch (fTrapType)
1018 {
1019 case 0: // General case
1020 {
1021 G4double dz = std::abs(p.z())-fDz;
1022 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1023 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1024 G4double dy = std::max(dz,std::max(dy1,dy2));
1025
1026 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1027 + fPlanes[2].c*p.z()+fPlanes[2].d;
1028 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1029 + fPlanes[3].c*p.z()+fPlanes[3].d;
1030 G4double dist = std::max(dy,std::max(dx1,dx2));
1031 return (dist < 0) ? -dist : 0.;
1032 }
1033 case 1: // YZ section is a rectangle
1034 {
1035 G4double dz = std::abs(p.z())-fDz;
1036 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1037 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1038 + fPlanes[2].c*p.z()+fPlanes[2].d;
1039 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1040 + fPlanes[3].c*p.z()+fPlanes[3].d;
1041 G4double dist = std::max(dy,std::max(dx1,dx2));
1042 return (dist < 0) ? -dist : 0.;
1043 }
1044 case 2: // YZ section is a rectangle and
1045 { // XZ section is an isosceles trapezoid
1046 G4double dz = std::abs(p.z())-fDz;
1047 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1048 G4double dx = fPlanes[3].a*std::abs(p.x())
1049 + fPlanes[3].c*p.z()+fPlanes[3].d;
1050 G4double dist = std::max(dy,dx);
1051 return (dist < 0) ? -dist : 0.;
1052 }
1053 case 3: // YZ section is a rectangle and
1054 { // XY section is an isosceles trapezoid
1055 G4double dz = std::abs(p.z())-fDz;
1056 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1057 G4double dx = fPlanes[3].a*std::abs(p.x())
1058 + fPlanes[3].b*p.y()+fPlanes[3].d;
1059 G4double dist = std::max(dy,dx);
1060 return (dist < 0) ? -dist : 0.;
1061 }
1062 }
1063 return 0.;
1064}
1065
1066//////////////////////////////////////////////////////////////////////////
1067//
1068// GetEntityType
1069
1070G4GeometryType G4Trap::GetEntityType() const
1071{
1072 return {"G4Trap"};
1073}
1074
1075//////////////////////////////////////////////////////////////////////////
1076//
1077// IsFaceted
1078
1079G4bool G4Trap::IsFaceted() const
1080{
1081 return true;
1082}
1083
1084//////////////////////////////////////////////////////////////////////////
1085//
1086// Make a clone of the object
1087//
1088G4VSolid* G4Trap::Clone() const
1089{
1090 return new G4Trap(*this);
1091}
1092
1093//////////////////////////////////////////////////////////////////////////
1094//
1095// Stream object contents to an output stream
1096
1097std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1098{
1099 G4double phi = GetPhi();
1100 G4double theta = GetTheta();
1101 G4double alpha1 = GetAlpha1();
1102 G4double alpha2 = GetAlpha2();
1103
1104 G4long oldprc = os.precision(16);
1105 os << "-----------------------------------------------------------\n"
1106 << " *** Dump for solid: " << GetName() << " ***\n"
1107 << " ===================================================\n"
1108 << " Solid type: G4Trap\n"
1109 << " Parameters:\n"
1110 << " half length Z: " << fDz/mm << " mm\n"
1111 << " half length Y, face -Dz: " << fDy1/mm << " mm\n"
1112 << " half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n"
1113 << " half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n"
1114 << " half length Y, face +Dz: " << fDy2/mm << " mm\n"
1115 << " half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n"
1116 << " half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n"
1117 << " theta: " << theta/degree << " degrees\n"
1118 << " phi: " << phi/degree << " degrees\n"
1119 << " alpha, face -Dz: " << alpha1/degree << " degrees\n"
1120 << " alpha, face +Dz: " << alpha2/degree << " degrees\n"
1121 << "-----------------------------------------------------------\n";
1122 os.precision(oldprc);
1123
1124 return os;
1125}
1126
1127//////////////////////////////////////////////////////////////////////////
1128//
1129// Compute vertices from planes
1130
1131void G4Trap::GetVertices(G4ThreeVector pt[8]) const
1132{
1133 for (G4int i=0; i<8; ++i)
1134 {
1135 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
1136 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
1137 G4double z = (i < 4) ? -fDz : fDz;
1138 G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b;
1139 G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z
1140 + fPlanes[ix].d)/fPlanes[ix].a;
1141 pt[i].set(x,y,z);
1142 }
1143}
1144
1145//////////////////////////////////////////////////////////////////////////
1146//
1147// Generate random point on the surface
1148
1149G4ThreeVector G4Trap::GetPointOnSurface() const
1150{
1151 // Set indeces
1152 constexpr G4int iface [6][4] =
1153 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
1154
1155 // Set vertices
1156 G4ThreeVector pt[8];
1157 GetVertices(pt);
1158
1159 // Select face
1160 //
1161 G4double select = fAreas[5]*G4QuickRand();
1162 G4int k = 5;
1163 k -= (G4int)(select <= fAreas[4]);
1164 k -= (G4int)(select <= fAreas[3]);
1165 k -= (G4int)(select <= fAreas[2]);
1166 k -= (G4int)(select <= fAreas[1]);
1167 k -= (G4int)(select <= fAreas[0]);
1168
1169 // Select sub-triangle
1170 //
1171 G4int i0 = iface[k][0];
1172 G4int i1 = iface[k][1];
1173 G4int i2 = iface[k][2];
1174 G4int i3 = iface[k][3];
1175 G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
1176 if (select > fAreas[k] - s2) { i0 = i2; }
1177
1178 // Generate point
1179 //
1180 G4double u = G4QuickRand();
1181 G4double v = G4QuickRand();
1182 if (u + v > 1.) { u = 1. - u; v = 1. - v; }
1183 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
1184}
1185
1186//////////////////////////////////////////////////////////////////////////
1187//
1188// Get volume
1189
1190G4double G4Trap::GetCubicVolume()
1191{
1192 if (fCubicVolume == 0)
1193 {
1194 G4AutoLock l(&trapMutex);
1195 G4ThreeVector pt[8];
1196 GetVertices(pt);
1197
1198 G4double dz = pt[4].z() - pt[0].z();
1199 G4double dy1 = pt[2].y() - pt[0].y();
1200 G4double dx1 = pt[1].x() - pt[0].x();
1201 G4double dx2 = pt[3].x() - pt[2].x();
1202 G4double dy2 = pt[6].y() - pt[4].y();
1203 G4double dx3 = pt[5].x() - pt[4].x();
1204 G4double dx4 = pt[7].x() - pt[6].x();
1205
1206 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) +
1207 (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125;
1208 l.unlock();
1209 }
1210 return fCubicVolume;
1211}
1212
1213//////////////////////////////////////////////////////////////////////////
1214//
1215// Get surface area
1216
1217G4double G4Trap::GetSurfaceArea()
1218{
1219 if (fSurfaceArea == 0)
1220 {
1221 G4AutoLock l(&trapMutex);
1222 G4ThreeVector pt[8];
1223 G4int iface [6][4] =
1224 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
1225
1226 GetVertices(pt);
1227 for (const auto & i : iface)
1228 {
1229 fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[i[0]],
1230 pt[i[1]],
1231 pt[i[2]],
1232 pt[i[3]]).mag();
1233 }
1234 l.unlock();
1235 }
1236 return fSurfaceArea;
1237}
1238
1239//////////////////////////////////////////////////////////////////////////
1240//
1241// Methods for visualisation
1242
1243void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1244{
1245 scene.AddSolid (*this);
1246}
1247
1248G4Polyhedron* G4Trap::CreatePolyhedron () const
1249{
1250 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1251 G4double alpha1 = std::atan(fTalpha1);
1252 G4double alpha2 = std::atan(fTalpha2);
1253 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1254 +fTthetaSphi*fTthetaSphi));
1255
1256 return new G4PolyhedronTrap(fDz, theta, phi,
1257 fDy1, fDx1, fDx2, alpha1,
1258 fDy2, fDx3, fDx4, alpha2);
1259}
1260
1261#endif
G4AutoLock
G4TemplateAutoLock< G4Mutex > G4AutoLock
Definition G4AutoLock.hh:582
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition G4BoundingEnvelope.hh:54
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4QuickRand
G4double G4QuickRand(uint32_t seed=0)
Definition G4QuickRand.hh:34
G4MUTEX_INITIALIZER
#define G4MUTEX_INITIALIZER
Definition G4Threading.hh:85
G4Mutex
std::mutex G4Mutex
Definition G4Threading.hh:81
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4long
long G4long
Definition G4Types.hh:87
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4GeometryType
G4String G4GeometryType
Definition G4VSolid.hh:70
alpha2
const G4double alpha2
Definition G4WentzelOKandVIxSection.cc:70
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4AffineTransform
G4AffineTransform is a class for geometric affine transformations. It supports efficient arbitrary ro...
Definition G4AffineTransform.hh:69
G4BoundingEnvelope
G4BoundingEnvelope is a helper class to facilitate calculation of the extent of a solid within the li...
Definition G4BoundingEnvelope.hh:65
G4BoundingEnvelope::BoundingBoxVsVoxelLimits
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Definition G4BoundingEnvelope.cc:182
G4BoundingEnvelope::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Definition G4BoundingEnvelope.cc:267
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition G4CSGSolid.hh:93
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition G4CSGSolid.hh:92
G4CSGSolid::fRebuildPolyhedron
G4bool fRebuildPolyhedron
Definition G4CSGSolid.hh:94
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition G4CSGSolid.cc:49
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition G4CSGSolid.cc:89
G4GeomTools::QuadAreaNormal
static G4ThreeVector QuadAreaNormal(const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D)
Definition G4GeomTools.cc:618
G4GeomTools::TriangleAreaNormal
static G4ThreeVector TriangleAreaNormal(const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C)
Definition G4GeomTools.cc:607
G4Trap::GetCubicVolume
G4double GetCubicVolume() override
Definition G4Trap.cc:1190
G4Trap::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep) override
Definition G4Trap.cc:487
G4Trap::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const override
Definition G4Trap.cc:538
G4Trap::IsFaceted
G4bool IsFaceted() const override
Definition G4Trap.cc:1079
G4Trap::GetEntityType
G4GeometryType GetEntityType() const override
Definition G4Trap.cc:1070
G4Trap::SetAllParameters
void SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
Definition G4Trap.cc:280
G4Trap::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const override
Definition G4Trap.cc:1149
G4Trap::GetAlpha2
G4double GetAlpha2() const
G4Trap::GetAlpha1
G4double GetAlpha1() const
G4Trap::GetTheta
G4double GetTheta() const
G4Trap::GetPhi
G4double GetPhi() const
G4Trap::SetCachedValues
void SetCachedValues()
Definition G4Trap.cc:432
G4Trap::Clone
G4VSolid * Clone() const override
Definition G4Trap.cc:1088
G4Trap::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const override
Definition G4Trap.cc:783
G4Trap::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const override
Definition G4Trap.cc:650
G4Trap::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const override
Definition G4Trap.cc:1097
G4Trap::Inside
EInside Inside(const G4ThreeVector &p) const override
Definition G4Trap.cc:587
G4Trap::GetSurfaceArea
G4double GetSurfaceArea() override
Definition G4Trap.cc:1217
G4Trap::GetZHalfLength
G4double GetZHalfLength() const
G4Trap::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const override
Definition G4Trap.cc:1248
G4Trap::MakePlanes
void MakePlanes()
Definition G4Trap.cc:333
G4Trap::G4Trap
G4Trap(const G4String &pName, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
Definition G4Trap.cc:66
G4Trap::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
Definition G4Trap.cc:912
G4Trap::MakePlane
G4bool MakePlane(const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, TrapSidePlane &plane)
Definition G4Trap.cc:400
G4Trap::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const override
Definition G4Trap.cc:1243
G4Trap::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
Definition G4Trap.cc:498
G4Trap::operator=
G4Trap & operator=(const G4Trap &rhs)
Definition G4Trap.cc:253
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4VPVParameterisation
G4VPVParameterisation ia an abstract base class for Parameterisation, able to compute the transformat...
Definition G4VPVParameterisation.hh:72
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:123
G4VPhysicalVolume
G4VPhysicalVolume is an abstract base class for the representation of a positioned volume....
Definition G4VPhysicalVolume.hh:85
G4VSolid::GetName
G4String GetName() const
G4VSolid::G4VSolid
G4VSolid(const G4String &name)
Definition G4VSolid.cc:59
G4VSolid::DumpInfo
void DumpInfo() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition G4VSolid.hh:418
G4VoxelLimits
G4VoxelLimits represents limitation/restrictions of space, where restrictions are only made perpendic...
Definition G4VoxelLimits.hh:58
EAxis
EAxis
Definition geomdefs.hh:54
EInside
EInside
Definition geomdefs.hh:67
kInside
@ kInside
Definition geomdefs.hh:70
kOutside
@ kOutside
Definition geomdefs.hh:68
kSurface
@ kSurface
Definition geomdefs.hh:69
TrapSidePlane::b
G4double b
Definition G4Trap.hh:90
TrapSidePlane::c
G4double c
Definition G4Trap.hh:90
TrapSidePlane::d
G4double d
Definition G4Trap.hh:90
TrapSidePlane::a
G4double a
Definition G4Trap.hh:90
DBL_EPSILON
#define DBL_EPSILON
Definition templates.hh:66
DBL_MAX
#define DBL_MAX
Definition templates.hh:62