| Geant4 Cross Reference |
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
53 using namespace CLHEP;
54
55 //////////////////////////////////////////////////////////////////////////
56 //
57 // Constructor - check and set half-widths as well as angles:
58 // final check of coplanarity
59
60 G4Trap::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
86 G4Trap::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
145 G4Trap::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
163 G4Trap::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
181 G4Trap::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
205 G4Trap::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 //
219 G4Trap::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
232 G4Trap::~G4Trap() = default;
233
234 //////////////////////////////////////////////////////////////////////////
235 //
236 // Copy constructor
237
238 G4Trap::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
253 G4Trap& 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
280 void 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
313 void 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
333 void 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
359 void 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
400 G4bool 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
432 void 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 // Get volume
485
486 G4double G4Trap::GetCubicVolume()
487 {
488 if (fCubicVolume == 0)
489 {
490 G4ThreeVector pt[8];
491 GetVertices(pt);
492
493 G4double dz = pt[4].z() - pt[0].z();
494 G4double dy1 = pt[2].y() - pt[0].y();
495 G4double dx1 = pt[1].x() - pt[0].x();
496 G4double dx2 = pt[3].x() - pt[2].x();
497 G4double dy2 = pt[6].y() - pt[4].y();
498 G4double dx3 = pt[5].x() - pt[4].x();
499 G4double dx4 = pt[7].x() - pt[6].x();
500
501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) +
502 (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125;
503 }
504 return fCubicVolume;
505 }
506
507 //////////////////////////////////////////////////////////////////////////
508 //
509 // Get surface area
510
511 G4double G4Trap::GetSurfaceArea()
512 {
513 if (fSurfaceArea == 0)
514 {
515 G4ThreeVector pt[8];
516 G4int iface [6][4] =
517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
518
519 GetVertices(pt);
520 for (const auto & i : iface)
521 {
522 fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[i[0]],
523 pt[i[1]],
524 pt[i[2]],
525 pt[i[3]]).mag();
526 }
527 }
528 return fSurfaceArea;
529 }
530
531 //////////////////////////////////////////////////////////////////////////
532 //
533 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification.
535
536 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
537 const G4int n,
538 const G4VPhysicalVolume* pRep )
539 {
540 p->ComputeDimensions(*this,n,pRep);
541 }
542
543 //////////////////////////////////////////////////////////////////////////
544 //
545 // Get bounding box
546
547 void G4Trap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
548 {
549 G4ThreeVector pt[8];
550 GetVertices(pt);
551
552 G4double xmin = kInfinity, xmax = -kInfinity;
553 G4double ymin = kInfinity, ymax = -kInfinity;
554 for (const auto & i : pt)
555 {
556 G4double x = i.x();
557 if (x < xmin) xmin = x;
558 if (x > xmax) xmax = x;
559 G4double y = i.y();
560 if (y < ymin) ymin = y;
561 if (y > ymax) ymax = y;
562 }
563
564 G4double dz = GetZHalfLength();
565 pMin.set(xmin,ymin,-dz);
566 pMax.set(xmax,ymax, dz);
567
568 // Check correctness of the bounding box
569 //
570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
571 {
572 std::ostringstream message;
573 message << "Bad bounding box (min >= max) for solid: "
574 << GetName() << " !"
575 << "\npMin = " << pMin
576 << "\npMax = " << pMax;
577 G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001",
578 JustWarning, message);
579 DumpInfo();
580 }
581 }
582
583 //////////////////////////////////////////////////////////////////////////
584 //
585 // Calculate extent under transform and specified limit
586
587 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588 const G4VoxelLimits& pVoxelLimit,
589 const G4AffineTransform& pTransform,
590 G4double& pMin, G4double& pMax) const
591 {
592 G4ThreeVector bmin, bmax;
593 G4bool exist;
594
595 // Check bounding box (bbox)
596 //
597 BoundingLimits(bmin,bmax);
598 G4BoundingEnvelope bbox(bmin,bmax);
599 #ifdef G4BBOX_EXTENT
600 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
601 #endif
602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
603 {
604 return exist = pMin < pMax;
605 }
606
607 // Set bounding envelope (benv) and calculate extent
608 //
609 G4ThreeVector pt[8];
610 GetVertices(pt);
611
612 G4ThreeVectorList baseA(4), baseB(4);
613 baseA[0] = pt[0];
614 baseA[1] = pt[1];
615 baseA[2] = pt[3];
616 baseA[3] = pt[2];
617
618 baseB[0] = pt[4];
619 baseB[1] = pt[5];
620 baseB[2] = pt[7];
621 baseB[3] = pt[6];
622
623 std::vector<const G4ThreeVectorList *> polygons(2);
624 polygons[0] = &baseA;
625 polygons[1] = &baseB;
626
627 G4BoundingEnvelope benv(bmin,bmax,polygons);
628 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
629 return exist;
630 }
631
632 //////////////////////////////////////////////////////////////////////////
633 //
634 // Return whether point is inside/outside/on_surface
635
636 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 {
638 switch (fTrapType)
639 {
640 case 0: // General case
641 {
642 G4double dz = std::abs(p.z())-fDz;
643 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
644 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
645 G4double dy = std::max(dz,std::max(dy1,dy2));
646
647 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
648 + fPlanes[2].c*p.z()+fPlanes[2].d;
649 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
650 + fPlanes[3].c*p.z()+fPlanes[3].d;
651 G4double dist = std::max(dy,std::max(dx1,dx2));
652
653 return (dist > halfCarTolerance) ? kOutside :
654 ((dist > -halfCarTolerance) ? kSurface : kInside);
655 }
656 case 1: // YZ section is a rectangle
657 {
658 G4double dz = std::abs(p.z())-fDz;
659 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
660 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
661 + fPlanes[2].c*p.z()+fPlanes[2].d;
662 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
663 + fPlanes[3].c*p.z()+fPlanes[3].d;
664 G4double dist = std::max(dy,std::max(dx1,dx2));
665
666 return (dist > halfCarTolerance) ? kOutside :
667 ((dist > -halfCarTolerance) ? kSurface : kInside);
668 }
669 case 2: // YZ section is a rectangle and
670 { // XZ section is an isosceles trapezoid
671 G4double dz = std::abs(p.z())-fDz;
672 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
673 G4double dx = fPlanes[3].a*std::abs(p.x())
674 + fPlanes[3].c*p.z()+fPlanes[3].d;
675 G4double dist = std::max(dy,dx);
676
677 return (dist > halfCarTolerance) ? kOutside :
678 ((dist > -halfCarTolerance) ? kSurface : kInside);
679 }
680 case 3: // YZ section is a rectangle and
681 { // XY section is an isosceles trapezoid
682 G4double dz = std::abs(p.z())-fDz;
683 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
684 G4double dx = fPlanes[3].a*std::abs(p.x())
685 + fPlanes[3].b*p.y()+fPlanes[3].d;
686 G4double dist = std::max(dy,dx);
687
688 return (dist > halfCarTolerance) ? kOutside :
689 ((dist > -halfCarTolerance) ? kSurface : kInside);
690 }
691 }
692 return kOutside;
693 }
694
695 //////////////////////////////////////////////////////////////////////////
696 //
697 // Determine side, and return corresponding normal
698
699 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 {
701 G4double nx = 0, ny = 0, nz = 0;
702 G4double dz = std::abs(p.z()) - fDz;
703 nz = std::copysign(G4double(std::abs(dz) <= halfCarTolerance), p.z());
704
705 switch (fTrapType)
706 {
707 case 0: // General case
708 {
709 for (G4int i=0; i<2; ++i)
710 {
711 G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
712 if (std::abs(dy) > halfCarTolerance) continue;
713 ny = fPlanes[i].b;
714 nz += fPlanes[i].c;
715 break;
716 }
717 for (G4int i=2; i<4; ++i)
718 {
719 G4double dx = fPlanes[i].a*p.x() +
720 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
721 if (std::abs(dx) > halfCarTolerance) continue;
722 nx = fPlanes[i].a;
723 ny += fPlanes[i].b;
724 nz += fPlanes[i].c;
725 break;
726 }
727 break;
728 }
729 case 1: // YZ section - rectangle
730 {
731 G4double dy = std::abs(p.y()) + fPlanes[1].d;
732 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
733 for (G4int i=2; i<4; ++i)
734 {
735 G4double dx = fPlanes[i].a*p.x() +
736 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
737 if (std::abs(dx) > halfCarTolerance) continue;
738 nx = fPlanes[i].a;
739 ny += fPlanes[i].b;
740 nz += fPlanes[i].c;
741 break;
742 }
743 break;
744 }
745 case 2: // YZ section - rectangle, XZ section - isosceles trapezoid
746 {
747 G4double dy = std::abs(p.y()) + fPlanes[1].d;
748 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
749 G4double dx = fPlanes[3].a*std::abs(p.x()) +
750 fPlanes[3].c*p.z() + fPlanes[3].d;
751 G4double k = std::abs(dx) <= halfCarTolerance;
752 nx = std::copysign(k, p.x())*fPlanes[3].a;
753 nz += k*fPlanes[3].c;
754 break;
755 }
756 case 3: // YZ section - rectangle, XY section - isosceles trapezoid
757 {
758 G4double dy = std::abs(p.y()) + fPlanes[1].d;
759 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
760 G4double dx = fPlanes[3].a*std::abs(p.x()) +
761 fPlanes[3].b*p.y() + fPlanes[3].d;
762 G4double k = std::abs(dx) <= halfCarTolerance;
763 nx = std::copysign(k, p.x())*fPlanes[3].a;
764 ny += k*fPlanes[3].b;
765 break;
766 }
767 }
768
769 // Return normal
770 //
771 G4double mag2 = nx*nx + ny*ny + nz*nz;
772 if (mag2 == 1) return { nx,ny,nz };
773 else if (mag2 != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
774 else
775 {
776 // Point is not on the surface
777 //
778 #ifdef G4CSGDEBUG
779 std::ostringstream message;
780 G4long oldprc = message.precision(16);
781 message << "Point p is not on surface (!?) of solid: "
782 << GetName() << G4endl;
783 message << "Position:\n";
784 message << " p.x() = " << p.x()/mm << " mm\n";
785 message << " p.y() = " << p.y()/mm << " mm\n";
786 message << " p.z() = " << p.z()/mm << " mm";
787 G4cout.precision(oldprc) ;
788 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
789 JustWarning, message );
790 DumpInfo();
791 #endif
792 return ApproxSurfaceNormal(p);
793 }
794 }
795
796 //////////////////////////////////////////////////////////////////////////
797 //
798 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface
800
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 {
803 G4double dist = -DBL_MAX;
804 G4int iside = 0;
805 for (G4int i=0; i<4; ++i)
806 {
807 G4double d = fPlanes[i].a*p.x() +
808 fPlanes[i].b*p.y() +
809 fPlanes[i].c*p.z() + fPlanes[i].d;
810 if (d > dist) { dist = d; iside = i; }
811 }
812
813 G4double distz = std::abs(p.z()) - fDz;
814 if (dist > distz)
815 return { fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c };
816 else
817 return { 0, 0, (G4double)((p.z() < 0) ? -1 : 1) };
818 }
819
820 //////////////////////////////////////////////////////////////////////////
821 //
822 // Calculate distance to shape from outside
823 // - return kInfinity if no intersection
824
825 G4double G4Trap::DistanceToIn(const G4ThreeVector& p,
826 const G4ThreeVector& v ) const
827 {
828 // Z intersections
829 //
830 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
831 return kInfinity;
832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
833 G4double dz = (invz < 0) ? fDz : -fDz;
834 G4double tzmin = (p.z() + dz)*invz;
835 G4double tzmax = (p.z() - dz)*invz;
836
837 // Y intersections
838 //
839 G4double tymin = 0, tymax = DBL_MAX;
840 G4int i = 0;
841 for ( ; i<2; ++i)
842 {
843 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
844 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
845 if (dist >= -halfCarTolerance)
846 {
847 if (cosa >= 0) return kInfinity;
848 G4double tmp = -dist/cosa;
849 if (tymin < tmp) tymin = tmp;
850 }
851 else if (cosa > 0)
852 {
853 G4double tmp = -dist/cosa;
854 if (tymax > tmp) tymax = tmp;
855 }
856 }
857
858 // Z intersections
859 //
860 G4double txmin = 0, txmax = DBL_MAX;
861 for ( ; i<4; ++i)
862 {
863 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
864 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() +
865 fPlanes[i].d;
866 if (dist >= -halfCarTolerance)
867 {
868 if (cosa >= 0) return kInfinity;
869 G4double tmp = -dist/cosa;
870 if (txmin < tmp) txmin = tmp;
871 }
872 else if (cosa > 0)
873 {
874 G4double tmp = -dist/cosa;
875 if (txmax > tmp) txmax = tmp;
876 }
877 }
878
879 // Find distance
880 //
881 G4double tmin = std::max(std::max(txmin,tymin),tzmin);
882 G4double tmax = std::min(std::min(txmax,tymax),tzmax);
883
884 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
885 return (tmin < halfCarTolerance ) ? 0. : tmin;
886 }
887
888 //////////////////////////////////////////////////////////////////////////
889 //
890 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside
893
894 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 {
896 switch (fTrapType)
897 {
898 case 0: // General case
899 {
900 G4double dz = std::abs(p.z())-fDz;
901 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
902 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
903 G4double dy = std::max(dz,std::max(dy1,dy2));
904
905 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
906 + fPlanes[2].c*p.z()+fPlanes[2].d;
907 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
908 + fPlanes[3].c*p.z()+fPlanes[3].d;
909 G4double dist = std::max(dy,std::max(dx1,dx2));
910 return (dist > 0) ? dist : 0.;
911 }
912 case 1: // YZ section is a rectangle
913 {
914 G4double dz = std::abs(p.z())-fDz;
915 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
916 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
917 + fPlanes[2].c*p.z()+fPlanes[2].d;
918 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
919 + fPlanes[3].c*p.z()+fPlanes[3].d;
920 G4double dist = std::max(dy,std::max(dx1,dx2));
921 return (dist > 0) ? dist : 0.;
922 }
923 case 2: // YZ section is a rectangle and
924 { // XZ section is an isosceles trapezoid
925 G4double dz = std::abs(p.z())-fDz;
926 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
927 G4double dx = fPlanes[3].a*std::abs(p.x())
928 + fPlanes[3].c*p.z()+fPlanes[3].d;
929 G4double dist = std::max(dy,dx);
930 return (dist > 0) ? dist : 0.;
931 }
932 case 3: // YZ section is a rectangle and
933 { // XY section is an isosceles trapezoid
934 G4double dz = std::abs(p.z())-fDz;
935 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
936 G4double dx = fPlanes[3].a*std::abs(p.x())
937 + fPlanes[3].b*p.y()+fPlanes[3].d;
938 G4double dist = std::max(dy,dx);
939 return (dist > 0) ? dist : 0.;
940 }
941 }
942 return 0.;
943 }
944
945 //////////////////////////////////////////////////////////////////////////
946 //
947 // Calculate distance to surface of shape from inside and
948 // find normal at exit point, if required
949 // - when leaving the surface, return 0
950
951 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952 const G4bool calcNorm,
953 G4bool* validNorm, G4ThreeVector* n) const
954 {
955 // Z intersections
956 //
957 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
958 {
959 if (calcNorm)
960 {
961 *validNorm = true;
962 n->set(0, 0, (p.z() < 0) ? -1 : 1);
963 }
964 return 0;
965 }
966 G4double vz = v.z();
967 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
968 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
969
970 // Y intersections
971 //
972 G4int i = 0;
973 for ( ; i<2; ++i)
974 {
975 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
976 if (cosa > 0)
977 {
978 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
979 if (dist >= -halfCarTolerance)
980 {
981 if (calcNorm)
982 {
983 *validNorm = true;
984 n->set(0, fPlanes[i].b, fPlanes[i].c);
985 }
986 return 0;
987 }
988 G4double tmp = -dist/cosa;
989 if (tmax > tmp) { tmax = tmp; iside = i; }
990 }
991 }
992
993 // X intersections
994 //
995 for ( ; i<4; ++i)
996 {
997 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
998 if (cosa > 0)
999 {
1000 G4double dist = fPlanes[i].a*p.x() +
1001 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
1002 if (dist >= -halfCarTolerance)
1003 {
1004 if (calcNorm)
1005 {
1006 *validNorm = true;
1007 n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
1008 }
1009 return 0;
1010 }
1011 G4double tmp = -dist/cosa;
1012 if (tmax > tmp) { tmax = tmp; iside = i; }
1013 }
1014 }
1015
1016 // Set normal, if required, and return distance
1017 //
1018 if (calcNorm)
1019 {
1020 *validNorm = true;
1021 if (iside < 0)
1022 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
1023 else
1024 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
1025 }
1026 return tmax;
1027 }
1028
1029 //////////////////////////////////////////////////////////////////////////
1030 //
1031 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside
1033
1034 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 {
1036 #ifdef G4CSGDEBUG
1037 if( Inside(p) == kOutside )
1038 {
1039 std::ostringstream message;
1040 G4long oldprc = message.precision(16);
1041 message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
1042 message << "Position:\n";
1043 message << " p.x() = " << p.x()/mm << " mm\n";
1044 message << " p.y() = " << p.y()/mm << " mm\n";
1045 message << " p.z() = " << p.z()/mm << " mm";
1046 G4cout.precision(oldprc);
1047 G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
1048 JustWarning, message );
1049 DumpInfo();
1050 }
1051 #endif
1052 switch (fTrapType)
1053 {
1054 case 0: // General case
1055 {
1056 G4double dz = std::abs(p.z())-fDz;
1057 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1058 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1059 G4double dy = std::max(dz,std::max(dy1,dy2));
1060
1061 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1062 + fPlanes[2].c*p.z()+fPlanes[2].d;
1063 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1064 + fPlanes[3].c*p.z()+fPlanes[3].d;
1065 G4double dist = std::max(dy,std::max(dx1,dx2));
1066 return (dist < 0) ? -dist : 0.;
1067 }
1068 case 1: // YZ section is a rectangle
1069 {
1070 G4double dz = std::abs(p.z())-fDz;
1071 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1072 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1073 + fPlanes[2].c*p.z()+fPlanes[2].d;
1074 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1075 + fPlanes[3].c*p.z()+fPlanes[3].d;
1076 G4double dist = std::max(dy,std::max(dx1,dx2));
1077 return (dist < 0) ? -dist : 0.;
1078 }
1079 case 2: // YZ section is a rectangle and
1080 { // XZ section is an isosceles trapezoid
1081 G4double dz = std::abs(p.z())-fDz;
1082 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1083 G4double dx = fPlanes[3].a*std::abs(p.x())
1084 + fPlanes[3].c*p.z()+fPlanes[3].d;
1085 G4double dist = std::max(dy,dx);
1086 return (dist < 0) ? -dist : 0.;
1087 }
1088 case 3: // YZ section is a rectangle and
1089 { // XY section is an isosceles trapezoid
1090 G4double dz = std::abs(p.z())-fDz;
1091 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1092 G4double dx = fPlanes[3].a*std::abs(p.x())
1093 + fPlanes[3].b*p.y()+fPlanes[3].d;
1094 G4double dist = std::max(dy,dx);
1095 return (dist < 0) ? -dist : 0.;
1096 }
1097 }
1098 return 0.;
1099 }
1100
1101 //////////////////////////////////////////////////////////////////////////
1102 //
1103 // GetEntityType
1104
1105 G4GeometryType G4Trap::GetEntityType() const
1106 {
1107 return {"G4Trap"};
1108 }
1109
1110 //////////////////////////////////////////////////////////////////////////
1111 //
1112 // IsFaceted
1113
1114 G4bool G4Trap::IsFaceted() const
1115 {
1116 return true;
1117 }
1118
1119 //////////////////////////////////////////////////////////////////////////
1120 //
1121 // Make a clone of the object
1122 //
1123 G4VSolid* G4Trap::Clone() const
1124 {
1125 return new G4Trap(*this);
1126 }
1127
1128 //////////////////////////////////////////////////////////////////////////
1129 //
1130 // Stream object contents to an output stream
1131
1132 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 {
1134 G4double phi = GetPhi();
1135 G4double theta = GetTheta();
1136 G4double alpha1 = GetAlpha1();
1137 G4double alpha2 = GetAlpha2();
1138
1139 G4long oldprc = os.precision(16);
1140 os << "-----------------------------------------------------------\n"
1141 << " *** Dump for solid: " << GetName() << " ***\n"
1142 << " ===================================================\n"
1143 << " Solid type: G4Trap\n"
1144 << " Parameters:\n"
1145 << " half length Z: " << fDz/mm << " mm\n"
1146 << " half length Y, face -Dz: " << fDy1/mm << " mm\n"
1147 << " half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n"
1148 << " half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n"
1149 << " half length Y, face +Dz: " << fDy2/mm << " mm\n"
1150 << " half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n"
1151 << " half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n"
1152 << " theta: " << theta/degree << " degrees\n"
1153 << " phi: " << phi/degree << " degrees\n"
1154 << " alpha, face -Dz: " << alpha1/degree << " degrees\n"
1155 << " alpha, face +Dz: " << alpha2/degree << " degrees\n"
1156 << "-----------------------------------------------------------\n";
1157 os.precision(oldprc);
1158
1159 return os;
1160 }
1161
1162 //////////////////////////////////////////////////////////////////////////
1163 //
1164 // Compute vertices from planes
1165
1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) const
1167 {
1168 for (G4int i=0; i<8; ++i)
1169 {
1170 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
1171 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
1172 G4double z = (i < 4) ? -fDz : fDz;
1173 G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b;
1174 G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z
1175 + fPlanes[ix].d)/fPlanes[ix].a;
1176 pt[i].set(x,y,z);
1177 }
1178 }
1179
1180 //////////////////////////////////////////////////////////////////////////
1181 //
1182 // Generate random point on the surface
1183
1184 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 {
1186 // Set indeces
1187 constexpr G4int iface [6][4] =
1188 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
1189
1190 // Set vertices
1191 G4ThreeVector pt[8];
1192 GetVertices(pt);
1193
1194 // Select face
1195 //
1196 G4double select = fAreas[5]*G4QuickRand();
1197 G4int k = 5;
1198 k -= (G4int)(select <= fAreas[4]);
1199 k -= (G4int)(select <= fAreas[3]);
1200 k -= (G4int)(select <= fAreas[2]);
1201 k -= (G4int)(select <= fAreas[1]);
1202 k -= (G4int)(select <= fAreas[0]);
1203
1204 // Select sub-triangle
1205 //
1206 G4int i0 = iface[k][0];
1207 G4int i1 = iface[k][1];
1208 G4int i2 = iface[k][2];
1209 G4int i3 = iface[k][3];
1210 G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
1211 if (select > fAreas[k] - s2) i0 = i2;
1212
1213 // Generate point
1214 //
1215 G4double u = G4QuickRand();
1216 G4double v = G4QuickRand();
1217 if (u + v > 1.) { u = 1. - u; v = 1. - v; }
1218 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
1219 }
1220
1221 //////////////////////////////////////////////////////////////////////////
1222 //
1223 // Methods for visualisation
1224
1225 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 {
1227 scene.AddSolid (*this);
1228 }
1229
1230 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 {
1232 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233 G4double alpha1 = std::atan(fTalpha1);
1234 G4double alpha2 = std::atan(fTalpha2);
1235 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1236 +fTthetaSphi*fTthetaSphi));
1237
1238 return new G4PolyhedronTrap(fDz, theta, phi,
1239 fDy1, fDx1, fDx2, alpha1,
1240 fDy2, fDx3, fDx4, alpha2);
1241 }
1242
1243 #endif
1244