ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Trd.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Trd.cc
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 G4Trd class
27 //
28 // 12.01.95 P.Kent: First version
29 // 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal
30 // 25.05.17 E.Tcherniaev: complete revision, speed-up
31 // --------------------------------------------------------------------
32 
33 #include "G4Trd.hh"
34 
35 #if !defined(G4GEOM_USE_UTRD)
36 
37 #include "G4GeomTools.hh"
38 
39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42 #include "Randomize.hh"
43 
44 #include "G4VPVParameterisation.hh"
45 
46 #include "G4VGraphicsScene.hh"
47 
48 using namespace CLHEP;
49 
51 //
52 // Constructor - set & check half widths
53 
54 G4Trd::G4Trd(const G4String& pName,
55  G4double pdx1, G4double pdx2,
56  G4double pdy1, G4double pdy2,
57  G4double pdz)
58  : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance),
59  fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(pdy2), fDz(pdz)
60 {
62  MakePlanes();
63 }
64 
66 //
67 // Fake default constructor - sets only member data and allocates memory
68 // for usage restricted to object persistency
69 //
70 G4Trd::G4Trd( __void__& a )
71  : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance),
72  fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fDz(1.)
73 {
74  MakePlanes();
75 }
76 
78 //
79 // Destructor
80 
82 {
83 }
84 
86 //
87 // Copy constructor
88 
89 G4Trd::G4Trd(const G4Trd& rhs)
90  : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
91  fDx1(rhs.fDx1), fDx2(rhs.fDx2),
92  fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz)
93 {
94  for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
95 }
96 
98 //
99 // Assignment operator
100 
102 {
103  // Check assignment to self
104  //
105  if (this == &rhs) { return *this; }
106 
107  // Copy base class data
108  //
110 
111  // Copy data
112  //
114  fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;
115  fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;
116  fDz = rhs.fDz;
117  for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
118 
119  return *this;
120 }
121 
123 //
124 // Set all parameters, as for constructor - set and check half-widths
125 
127  G4double pdy1, G4double pdy2, G4double pdz)
128 {
129  // Reset data of the base class
130  fCubicVolume = 0.;
131  fSurfaceArea = 0.;
132  fRebuildPolyhedron = true;
133 
134  // Set parameters
135  fDx1 = pdx1; fDx2 = pdx2;
136  fDy1 = pdy1; fDy2 = pdy2;
137  fDz = pdz;
138 
139  CheckParameters();
140  MakePlanes();
141 }
142 
144 //
145 // Check dimensions
146 
148 {
149  G4double dmin = 2*kCarTolerance;
150  if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy2 < 0 || fDz < dmin) ||
151  (fDx1 < dmin && fDx2 < dmin) ||
152  (fDy1 < dmin && fDy2 < dmin))
153  {
154  std::ostringstream message;
155  message << "Invalid (too small or negative) dimensions for Solid: "
156  << GetName()
157  << "\n X - " << fDx1 << ", " << fDx2
158  << "\n Y - " << fDy1 << ", " << fDy2
159  << "\n Z - " << fDz;
160  G4Exception("G4Trd::CheckParameters()", "GeomSolids0002",
161  FatalException, message);
162  }
163 }
164 
166 //
167 // Set side planes
168 
170 {
171  G4double dx = fDx1 - fDx2;
172  G4double dy = fDy1 - fDy2;
173  G4double dz = 2*fDz;
174  G4double magx = std::sqrt(dx*dx + dz*dz);
175  G4double magy = std::sqrt(dy*dy + dz*dz);
176 
177  // Set -Y & +Y planes
178  //
179  fPlanes[0].a = 0.;
180  fPlanes[0].b = -dz/magy;
181  fPlanes[0].c = dy/magy;
182  fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0].c*fDz;
183 
184  fPlanes[1].a = fPlanes[0].a;
185  fPlanes[1].b = -fPlanes[0].b;
186  fPlanes[1].c = fPlanes[0].c;
187  fPlanes[1].d = fPlanes[0].d;
188 
189  // Set -X & +X planes
190  //
191  fPlanes[2].a = -dz/magx;
192  fPlanes[2].b = 0.;
193  fPlanes[2].c = dx/magx;
194  fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2].c*fDz;
195 
196  fPlanes[3].a = -fPlanes[2].a;
197  fPlanes[3].b = fPlanes[2].b;
198  fPlanes[3].c = fPlanes[2].c;
199  fPlanes[3].d = fPlanes[2].d;
200 }
201 
203 //
204 // Get volume
205 
207 {
208  if (fCubicVolume == 0.)
209  {
210  fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+fDy2) +
211  (fDx2-fDx1)*(fDy2-fDy1)/3 );
212  }
213  return fCubicVolume;
214 }
215 
217 //
218 // Get surface area
219 
221 {
222  if (fSurfaceArea == 0.)
223  {
224  fSurfaceArea =
225  4*(fDx1*fDy1+fDx2*fDy2) +
226  2*(fDy1+fDy2)*std::hypot(fDx1-fDx2,2*fDz) +
227  2*(fDx1+fDx2)*std::hypot(fDy1-fDy2,2*fDz);
228  }
229  return fSurfaceArea;
230 }
231 
233 //
234 // Dispatch to parameterisation for replication mechanism dimension
235 // computation & modification
236 
238  const G4int n,
239  const G4VPhysicalVolume* pRep )
240 {
241  p->ComputeDimensions(*this,n,pRep);
242 }
243 
245 //
246 // Get bounding box
247 
249 {
250  G4double dx1 = GetXHalfLength1();
251  G4double dx2 = GetXHalfLength2();
252  G4double dy1 = GetYHalfLength1();
253  G4double dy2 = GetYHalfLength2();
255 
256  G4double xmax = std::max(dx1,dx2);
257  G4double ymax = std::max(dy1,dy2);
258  pMin.set(-xmax,-ymax,-dz);
259  pMax.set( xmax, ymax, dz);
260 
261  // Check correctness of the bounding box
262  //
263  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
264  {
265  std::ostringstream message;
266  message << "Bad bounding box (min >= max) for solid: "
267  << GetName() << " !"
268  << "\npMin = " << pMin
269  << "\npMax = " << pMax;
270  G4Exception("G4Trd::BoundingLimits()", "GeomMgt0001", JustWarning, message);
271  DumpInfo();
272  }
273 }
274 
276 //
277 // Calculate extent under transform and specified limit
278 
280  const G4VoxelLimits& pVoxelLimit,
281  const G4AffineTransform& pTransform,
282  G4double& pMin, G4double& pMax ) const
283 {
284  G4ThreeVector bmin, bmax;
285  G4bool exist;
286 
287  // Check bounding box (bbox)
288  //
289  BoundingLimits(bmin,bmax);
290  G4BoundingEnvelope bbox(bmin,bmax);
291 #ifdef G4BBOX_EXTENT
292  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
293 #endif
294  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
295  {
296  return exist = (pMin < pMax) ? true : false;
297  }
298 
299  // Set bounding envelope (benv) and calculate extent
300  //
301  G4double dx1 = GetXHalfLength1();
302  G4double dx2 = GetXHalfLength2();
303  G4double dy1 = GetYHalfLength1();
304  G4double dy2 = GetYHalfLength2();
306 
307  G4ThreeVectorList baseA(4), baseB(4);
308  baseA[0].set(-dx1,-dy1,-dz);
309  baseA[1].set( dx1,-dy1,-dz);
310  baseA[2].set( dx1, dy1,-dz);
311  baseA[3].set(-dx1, dy1,-dz);
312  baseB[0].set(-dx2,-dy2, dz);
313  baseB[1].set( dx2,-dy2, dz);
314  baseB[2].set( dx2, dy2, dz);
315  baseB[3].set(-dx2, dy2, dz);
316 
317  std::vector<const G4ThreeVectorList *> polygons(2);
318  polygons[0] = &baseA;
319  polygons[1] = &baseB;
320 
321  G4BoundingEnvelope benv(bmin,bmax,polygons);
322  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
323  return exist;
324 }
325 
327 //
328 // Return whether point inside/outside/on surface, using tolerance
329 
331 {
332  G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
333  G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
334  G4double dxy = std::max(dx,dy);
335 
336  G4double dz = std::abs(p.z())-fDz;
337  G4double dist = std::max(dz,dxy);
338 
339  if (dist > halfCarTolerance) return kOutside;
340  return (dist > -halfCarTolerance) ? kSurface : kInside;
341 }
342 
344 //
345 // Determine side where point is, and return corresponding normal
346 
348 {
349  G4int nsurf = 0; // number of surfaces where p is placed
350 
351  // Check Z faces
352  //
353  G4double nz = 0;
354  G4double dz = std::abs(p.z()) - fDz;
355  if (std::abs(dz) <= halfCarTolerance)
356  {
357  nz = (p.z() < 0) ? -1 : 1;
358  ++nsurf;
359  }
360 
361  // Check Y faces
362  //
363  G4double ny = 0;
364  G4double dy1 = fPlanes[0].b*p.y();
365  G4double dy2 = fPlanes[0].c*p.z() + fPlanes[0].d;
366  if (std::abs(dy2 + dy1) <= halfCarTolerance)
367  {
368  ny += fPlanes[0].b;
369  nz += fPlanes[0].c;
370  ++nsurf;
371  }
372  if (std::abs(dy2 - dy1) <= halfCarTolerance)
373  {
374  ny += fPlanes[1].b;
375  nz += fPlanes[1].c;
376  ++nsurf;
377  }
378 
379  // Check X faces
380  //
381  G4double nx = 0;
382  G4double dx1 = fPlanes[2].a*p.x();
383  G4double dx2 = fPlanes[2].c*p.z() + fPlanes[2].d;
384  if (std::abs(dx2 + dx1) <= halfCarTolerance)
385  {
386  nx += fPlanes[2].a;
387  nz += fPlanes[2].c;
388  ++nsurf;
389  }
390  if (std::abs(dx2 - dx1) <= halfCarTolerance)
391  {
392  nx += fPlanes[3].a;
393  nz += fPlanes[3].c;
394  ++nsurf;
395  }
396 
397  // Return normal
398  //
399  if (nsurf == 1) return G4ThreeVector(nx,ny,nz);
400  else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
401  else
402  {
403  // Point is not on the surface
404  //
405 #ifdef G4CSGDEBUG
406  std::ostringstream message;
407  G4int oldprc = message.precision(16);
408  message << "Point p is not on surface (!?) of solid: "
409  << GetName() << G4endl;
410  message << "Position:\n";
411  message << " p.x() = " << p.x()/mm << " mm\n";
412  message << " p.y() = " << p.y()/mm << " mm\n";
413  message << " p.z() = " << p.z()/mm << " mm";
414  G4cout.precision(oldprc) ;
415  G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002",
416  JustWarning, message );
417  DumpInfo();
418 #endif
419  return ApproxSurfaceNormal(p);
420  }
421 }
422 
424 //
425 // Algorithm for SurfaceNormal() following the original specification
426 // for points not on the surface
427 
429 {
430  G4double dist = -DBL_MAX;
431  G4int iside = 0;
432  for (G4int i=0; i<4; ++i)
433  {
434  G4double d = fPlanes[i].a*p.x() +
435  fPlanes[i].b*p.y() +
436  fPlanes[i].c*p.z() + fPlanes[i].d;
437  if (d > dist) { dist = d; iside = i; }
438  }
439 
440  G4double distz = std::abs(p.z()) - fDz;
441  if (dist > distz)
442  return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
443  else
444  return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
445 }
446 
448 //
449 // Calculate distance to shape from outside
450 // - return kInfinity if no intersection
451 
453  const G4ThreeVector& v ) const
454 {
455  // Z intersections
456  //
457  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
458  return kInfinity;
459  G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
460  G4double dz = (invz < 0) ? fDz : -fDz;
461  G4double tzmin = (p.z() + dz)*invz;
462  G4double tzmax = (p.z() - dz)*invz;
463 
464  // Y intersections
465  //
466  G4double tmin0 = tzmin, tmax0 = tzmax;
467  G4double ya = fPlanes[0].b*v.y(), yb = fPlanes[0].c*v.z();
468  G4double yc = fPlanes[0].b*p.y(), yd = fPlanes[0].c*p.z()+fPlanes[0].d;
469  G4double cos0 = yb + ya;
470  G4double dis0 = yd + yc;
471  if (dis0 >= -halfCarTolerance)
472  {
473  if (cos0 >= 0) return kInfinity;
474  G4double tmp = -dis0/cos0;
475  if (tmin0 < tmp) tmin0 = tmp;
476  }
477  else if (cos0 > 0)
478  {
479  G4double tmp = -dis0/cos0;
480  if (tmax0 > tmp) tmax0 = tmp;
481  }
482 
483  G4double tmin1 = tmin0, tmax1 = tmax0;
484  G4double cos1 = yb - ya;
485  G4double dis1 = yd - yc;
486  if (dis1 >= -halfCarTolerance)
487  {
488  if (cos1 >= 0) return kInfinity;
489  G4double tmp = -dis1/cos1;
490  if (tmin1 < tmp) tmin1 = tmp;
491  }
492  else if (cos1 > 0)
493  {
494  G4double tmp = -dis1/cos1;
495  if (tmax1 > tmp) tmax1 = tmp;
496  }
497 
498  // X intersections
499  //
500  G4double tmin2 = tmin1, tmax2 = tmax1;
501  G4double xa = fPlanes[2].a*v.x(), xb = fPlanes[2].c*v.z();
502  G4double xc = fPlanes[2].a*p.x(), xd = fPlanes[2].c*p.z()+fPlanes[2].d;
503  G4double cos2 = xb + xa;
504  G4double dis2 = xd + xc;
505  if (dis2 >= -halfCarTolerance)
506  {
507  if (cos2 >= 0) return kInfinity;
508  G4double tmp = -dis2/cos2;
509  if (tmin2 < tmp) tmin2 = tmp;
510  }
511  else if (cos2 > 0)
512  {
513  G4double tmp = -dis2/cos2;
514  if (tmax2 > tmp) tmax2 = tmp;
515  }
516 
517  G4double tmin3 = tmin2, tmax3 = tmax2;
518  G4double cos3 = xb - xa;
519  G4double dis3 = xd - xc;
520  if (dis3 >= -halfCarTolerance)
521  {
522  if (cos3 >= 0) return kInfinity;
523  G4double tmp = -dis3/cos3;
524  if (tmin3 < tmp) tmin3 = tmp;
525  }
526  else if (cos3 > 0)
527  {
528  G4double tmp = -dis3/cos3;
529  if (tmax3 > tmp) tmax3 = tmp;
530  }
531 
532  // Find distance
533  //
534  G4double tmin = tmin3, tmax = tmax3;
535  if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
536  return (tmin < halfCarTolerance ) ? 0. : tmin;
537 }
538 
540 //
541 // Calculate exact shortest distance to any boundary from outside
542 // This is the best fast estimation of the shortest distance to trap
543 // - returns 0 if point is inside
544 
546 {
547  G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
548  G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
549  G4double dxy = std::max(dx,dy);
550 
551  G4double dz = std::abs(p.z())-fDz;
552  G4double dist = std::max(dz,dxy);
553 
554  return (dist > 0) ? dist : 0.;
555 }
556 
558 //
559 // Calculate distance to surface of shape from inside and
560 // find normal at exit point, if required
561 // - when leaving the surface, return 0
562 
564  const G4bool calcNorm,
565  G4bool* validNorm, G4ThreeVector* n) const
566 {
567  // Z intersections
568  //
569  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
570  {
571  if (calcNorm)
572  {
573  *validNorm = true;
574  n->set(0, 0, (p.z() < 0) ? -1 : 1);
575  }
576  return 0;
577  }
578  G4double vz = v.z();
579  G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
580  G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
581 
582  // Y intersections
583  //
584  G4int i = 0;
585  for ( ; i<2; ++i)
586  {
587  G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
588  if (cosa > 0)
589  {
590  G4double dist = fPlanes[i].b*p.y()+fPlanes[i].c*p.z()+fPlanes[i].d;
591  if (dist >= -halfCarTolerance)
592  {
593  if (calcNorm)
594  {
595  *validNorm = true;
596  n->set(0, fPlanes[i].b, fPlanes[i].c);
597  }
598  return 0;
599  }
600  G4double tmp = -dist/cosa;
601  if (tmax > tmp) { tmax = tmp; iside = i; }
602  }
603  }
604 
605  // X intersections
606  //
607  for ( ; i<4; ++i)
608  {
609  G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].c*v.z();
610  if (cosa > 0)
611  {
612  G4double dist = fPlanes[i].a*p.x()+fPlanes[i].c*p.z()+fPlanes[i].d;
613  if (dist >= -halfCarTolerance)
614  {
615  if (calcNorm)
616  {
617  *validNorm = true;
618  n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
619  }
620  return 0;
621  }
622  G4double tmp = -dist/cosa;
623  if (tmax > tmp) { tmax = tmp; iside = i; }
624  }
625  }
626 
627  // Set normal, if required, and return distance
628  //
629  if (calcNorm)
630  {
631  *validNorm = true;
632  if (iside < 0)
633  n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
634  else
635  n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
636  }
637  return tmax;
638 }
639 
641 //
642 // Calculate exact shortest distance to any boundary from inside
643 // - returns 0 if point is outside
644 
646 {
647 #ifdef G4CSGDEBUG
648  if( Inside(p) == kOutside )
649  {
650  std::ostringstream message;
651  G4int oldprc = message.precision(16);
652  message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
653  message << "Position:\n";
654  message << " p.x() = " << p.x()/mm << " mm\n";
655  message << " p.y() = " << p.y()/mm << " mm\n";
656  message << " p.z() = " << p.z()/mm << " mm";
657  G4cout.precision(oldprc) ;
658  G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002",
659  JustWarning, message );
660  DumpInfo();
661  }
662 #endif
663  G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
664  G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
665  G4double dxy = std::max(dx,dy);
666 
667  G4double dz = std::abs(p.z())-fDz;
668  G4double dist = std::max(dz,dxy);
669 
670  return (dist < 0) ? -dist : 0.;
671 }
672 
674 //
675 // GetEntityType
676 
678 {
679  return G4String("G4Trd");
680 }
681 
683 //
684 // Make a clone of the object
685 //
687 {
688  return new G4Trd(*this);
689 }
690 
692 //
693 // Stream object contents to an output stream
694 
695 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
696 {
697  G4int oldprc = os.precision(16);
698  os << "-----------------------------------------------------------\n"
699  << " *** Dump for solid - " << GetName() << " ***\n"
700  << " ===================================================\n"
701  << " Solid type: G4Trd\n"
702  << " Parameters: \n"
703  << " half length X, surface -dZ: " << fDx1/mm << " mm \n"
704  << " half length X, surface +dZ: " << fDx2/mm << " mm \n"
705  << " half length Y, surface -dZ: " << fDy1/mm << " mm \n"
706  << " half length Y, surface +dZ: " << fDy2/mm << " mm \n"
707  << " half length Z : " << fDz/mm << " mm \n"
708  << "-----------------------------------------------------------\n";
709  os.precision(oldprc);
710 
711  return os;
712 }
713 
715 //
716 // Return a point randomly and uniformly selected on the solid surface
717 
719 {
720  // Set vertices
721  //
722  G4ThreeVector pt[8];
723  pt[0].set(-fDx1,-fDy1,-fDz);
724  pt[1].set( fDx1,-fDy1,-fDz);
725  pt[2].set(-fDx1, fDy1,-fDz);
726  pt[3].set( fDx1, fDy1,-fDz);
727  pt[4].set(-fDx2,-fDy2, fDz);
728  pt[5].set( fDx2,-fDy2, fDz);
729  pt[6].set(-fDx2, fDy2, fDz);
730  pt[7].set( fDx2, fDy2, fDz);
731 
732  // Set faces (-Z, -Y, +Y, -X, +X, +Z)
733  //
734  G4int iface [6][4] =
735  { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
736 
737  // Set areas
738  //
739  G4double sxz = (fDy1 + fDy2)*std::hypot(fDx1 - fDx2, 2*fDz);
740  G4double syz = (fDx1 + fDx2)*std::hypot(fDy1 - fDy2, 2*fDz);
741  G4double sface[6] = { 4*fDx1*fDy1, syz, syz, sxz, sxz, 4*fDx2*fDy2 };
742  for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; }
743 
744  // Select face
745  //
746  G4double select = sface[5]*G4UniformRand();
747  G4int k = 5;
748  if (select <= sface[4]) k = 4;
749  if (select <= sface[3]) k = 3;
750  if (select <= sface[2]) k = 2;
751  if (select <= sface[1]) k = 1;
752  if (select <= sface[0]) k = 0;
753 
754  // Select sub-triangle
755  //
756  G4int i0 = iface[k][0];
757  G4int i1 = iface[k][1];
758  G4int i2 = iface[k][2];
759  G4int i3 = iface[k][3];
760  G4double s1 = G4GeomTools::TriangleAreaNormal(pt[i0],pt[i1],pt[i3]).mag();
761  G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
762  if ((s1+s2)*G4UniformRand() > s1) i0 = i2;
763 
764  // Generate point
765  //
768  if (u + v > 1.) { u = 1. - u; v = 1. - v; }
769  return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
770 }
771 
773 //
774 // Methods for visualisation
775 
777 {
778  scene.AddSolid (*this);
779 }
780 
782 {
783  return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
784 }
785 
786 #endif