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