ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Trap.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Trap.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 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 "Randomize.hh"
49 
50 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh"
52 
53 using namespace CLHEP;
54 
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 
77  MakePlanes();
78 }
79 
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 
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;
154 
155  CheckParameters();
156  MakePlanes();
157 }
158 
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 
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 
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 
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 
229 //
230 // Destructor
231 
233 {
234 }
235 
237 //
238 // Copy constructor
239 
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  fTrapType = rhs.fTrapType;
248 }
249 
251 //
252 // Assignment operator
253 
255 {
256  // Check assignment to self
257  //
258  if (this == &rhs) { return *this; }
259 
260  // Copy base class data
261  //
263 
264  // Copy data
265  //
268  fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
269  fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
270  for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
271  fTrapType = rhs.fTrapType;
272  return *this;
273 }
274 
276 //
277 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanarity
279 
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 
310 //
311 // Check length parameters
312 
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 
330 //
331 // Compute vertices and set side planes
332 
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 
356 //
357 // Set side planes, check planarity
358 
360 {
361  G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} };
362  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
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  // Define type of trapezoid
389  fTrapType = 0;
390  if (fPlanes[0].b == -1 && fPlanes[1].b == 1 &&
391  std::abs(fPlanes[0].a) < DBL_EPSILON &&
392  std::abs(fPlanes[0].c) < DBL_EPSILON &&
393  std::abs(fPlanes[1].a) < DBL_EPSILON &&
394  std::abs(fPlanes[1].c) < DBL_EPSILON)
395  {
396  fTrapType = 1; // YZ section is a rectangle ...
397  if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
398  std::abs(fPlanes[2].c - fPlanes[3].c) < DBL_EPSILON &&
399  fPlanes[2].b == 0 &&
400  fPlanes[3].b == 0)
401  {
402  fTrapType = 2; // ... and XZ section is a isosceles trapezoid
403  fPlanes[2].a = -fPlanes[3].a;
404  fPlanes[2].c = fPlanes[3].c;
405  }
406  if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
407  std::abs(fPlanes[2].b - fPlanes[3].b) < DBL_EPSILON &&
408  fPlanes[2].c == 0 &&
409  fPlanes[3].c == 0)
410  {
411  fTrapType = 3; // ... and XY section is a isosceles trapezoid
412  fPlanes[2].a = -fPlanes[3].a;
413  fPlanes[2].b = fPlanes[3].b;
414  }
415  }
416 }
417 
419 //
420 // Calculate the coef's of the plane p1->p2->p3->p4->p1
421 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed
422 // from infront of the plane (i.e. from normal direction).
423 //
424 // Return true if the points are coplanar, false otherwise
425 
427  const G4ThreeVector& p2,
428  const G4ThreeVector& p3,
429  const G4ThreeVector& p4,
431 {
432  G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit();
433  if (std::abs(normal.x()) < DBL_EPSILON) normal.setX(0);
434  if (std::abs(normal.y()) < DBL_EPSILON) normal.setY(0);
435  if (std::abs(normal.z()) < DBL_EPSILON) normal.setZ(0);
436  normal = normal.unit();
437 
438  G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25;
439  plane.a = normal.x();
440  plane.b = normal.y();
441  plane.c = normal.z();
442  plane.d = -normal.dot(centre);
443 
444  // compute distances and check planarity
445  G4double d1 = std::abs(normal.dot(p1) + plane.d);
446  G4double d2 = std::abs(normal.dot(p2) + plane.d);
447  G4double d3 = std::abs(normal.dot(p3) + plane.d);
448  G4double d4 = std::abs(normal.dot(p4) + plane.d);
449  G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4);
450 
451  return (dmax > 1000 * kCarTolerance) ? false : true;
452 }
453 
455 //
456 // Get volume
457 
459 {
460  if (fCubicVolume == 0)
461  {
462  G4ThreeVector pt[8];
463  GetVertices(pt);
464 
465  G4double dz = pt[4].z() - pt[0].z();
466  G4double dy1 = pt[2].y() - pt[0].y();
467  G4double dx1 = pt[1].x() - pt[0].x();
468  G4double dx2 = pt[3].x() - pt[2].x();
469  G4double dy2 = pt[6].y() - pt[4].y();
470  G4double dx3 = pt[5].x() - pt[4].x();
471  G4double dx4 = pt[7].x() - pt[6].x();
472 
473  fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) +
474  (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125;
475  }
476  return fCubicVolume;
477 }
478 
480 //
481 // Get surface area
482 
484 {
485  if (fSurfaceArea == 0)
486  {
487  G4ThreeVector pt[8];
488  G4int iface [6][4] =
489  { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
490 
491  GetVertices(pt);
492  for (G4int i=0; i<6; ++i)
493  {
494  fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
495  pt[iface[i][1]],
496  pt[iface[i][2]],
497  pt[iface[i][3]]).mag();
498  }
499  }
500  return fSurfaceArea;
501 }
502 
504 //
505 // Dispatch to parameterisation for replication mechanism dimension
506 // computation & modification.
507 
509  const G4int n,
510  const G4VPhysicalVolume* pRep )
511 {
512  p->ComputeDimensions(*this,n,pRep);
513 }
514 
516 //
517 // Get bounding box
518 
520 {
521  G4ThreeVector pt[8];
522  GetVertices(pt);
523 
526  for (G4int i=0; i<8; ++i)
527  {
528  G4double x = pt[i].x();
529  if (x < xmin) xmin = x;
530  if (x > xmax) xmax = x;
531  G4double y = pt[i].y();
532  if (y < ymin) ymin = y;
533  if (y > ymax) ymax = y;
534  }
535 
537  pMin.set(xmin,ymin,-dz);
538  pMax.set(xmax,ymax, dz);
539 
540  // Check correctness of the bounding box
541  //
542  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
543  {
544  std::ostringstream message;
545  message << "Bad bounding box (min >= max) for solid: "
546  << GetName() << " !"
547  << "\npMin = " << pMin
548  << "\npMax = " << pMax;
549  G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001",
550  JustWarning, message);
551  DumpInfo();
552  }
553 }
554 
556 //
557 // Calculate extent under transform and specified limit
558 
560  const G4VoxelLimits& pVoxelLimit,
561  const G4AffineTransform& pTransform,
562  G4double& pMin, G4double& pMax) const
563 {
564  G4ThreeVector bmin, bmax;
565  G4bool exist;
566 
567  // Check bounding box (bbox)
568  //
569  BoundingLimits(bmin,bmax);
570  G4BoundingEnvelope bbox(bmin,bmax);
571 #ifdef G4BBOX_EXTENT
572  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
573 #endif
574  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
575  {
576  return exist = (pMin < pMax) ? true : false;
577  }
578 
579  // Set bounding envelope (benv) and calculate extent
580  //
581  G4ThreeVector pt[8];
582  GetVertices(pt);
583 
584  G4ThreeVectorList baseA(4), baseB(4);
585  baseA[0] = pt[0];
586  baseA[1] = pt[1];
587  baseA[2] = pt[3];
588  baseA[3] = pt[2];
589 
590  baseB[0] = pt[4];
591  baseB[1] = pt[5];
592  baseB[2] = pt[7];
593  baseB[3] = pt[6];
594 
595  std::vector<const G4ThreeVectorList *> polygons(2);
596  polygons[0] = &baseA;
597  polygons[1] = &baseB;
598 
599  G4BoundingEnvelope benv(bmin,bmax,polygons);
600  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
601  return exist;
602 }
603 
605 //
606 // Return whether point is inside/outside/on_surface
607 
609 {
610  G4double dz = std::abs(p.z())-fDz;
611  if (dz > halfCarTolerance) return kOutside;
612 
613  switch (fTrapType)
614  {
615  case 0: // General case
616  {
617  G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
618  G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
619  G4double dy = std::max(dz,std::max(dy1,dy2));
620 
621  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
622  + fPlanes[2].c*p.z()+fPlanes[2].d;
623  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
624  + fPlanes[3].c*p.z()+fPlanes[3].d;
625  G4double dist = std::max(dy,std::max(dx1,dx2));
626 
627  if (dist > halfCarTolerance) return kOutside;
628  return (dist > -halfCarTolerance) ? kSurface : kInside;
629  }
630  case 1: // YZ section is a rectangle
631  {
632  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
633  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
634  + fPlanes[2].c*p.z()+fPlanes[2].d;
635  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
636  + fPlanes[3].c*p.z()+fPlanes[3].d;
637  G4double dist = std::max(dy,std::max(dx1,dx2));
638 
639  if (dist > halfCarTolerance) return kOutside;
640  return (dist > -halfCarTolerance) ? kSurface : kInside;
641  }
642  case 2: // YZ section is a rectangle and
643  { // XZ section is an isosceles trapezoid
644  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
645  G4double dx = fPlanes[3].a*std::abs(p.x())
646  + fPlanes[3].c*p.z()+fPlanes[3].d;
647  G4double dist = std::max(dy,dx);
648 
649  if (dist > halfCarTolerance) return kOutside;
650  return (dist > -halfCarTolerance) ? kSurface : kInside;
651  }
652  case 3: // YZ section is a rectangle and
653  { // XY section is an isosceles trapezoid
654  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
655  G4double dx = fPlanes[3].a*std::abs(p.x())
656  + fPlanes[3].b*p.y()+fPlanes[3].d;
657  G4double dist = std::max(dy,dx);
658 
659  if (dist > halfCarTolerance) return kOutside;
660  return (dist > -halfCarTolerance) ? kSurface : kInside;
661  }
662  }
663  return kOutside;
664 }
665 
667 //
668 // Determine side, and return corresponding normal
669 
671 {
672  G4int nsurf = 0; // number of surfaces where p is placed
673  G4double nx = 0, ny = 0, nz = 0;
674  G4double dz = std::abs(p.z()) - fDz;
675  if (std::abs(dz) <= halfCarTolerance)
676  {
677  nz = (p.z() < 0) ? -1 : 1;
678  ++nsurf;
679  }
680 
681  switch (fTrapType)
682  {
683  case 0: // General case
684  {
685  for (G4int i=0; i<2; ++i)
686  {
687  G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
688  if (std::abs(dy) > halfCarTolerance) continue;
689  ny = fPlanes[i].b;
690  nz += fPlanes[i].c;
691  ++nsurf;
692  break;
693  }
694  for (G4int i=2; i<4; ++i)
695  {
696  G4double dx = fPlanes[i].a*p.x() +
697  fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
698  if (std::abs(dx) > halfCarTolerance) continue;
699  nx = fPlanes[i].a;
700  ny += fPlanes[i].b;
701  nz += fPlanes[i].c;
702  ++nsurf;
703  break;
704  }
705  break;
706  }
707  case 1: // YZ section is a rectangle
708  {
709  G4double dy = std::abs(p.y()) + fPlanes[1].d;
710  if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
711  for (G4int i=2; i<4; ++i)
712  {
713  G4double dx = fPlanes[i].a*p.x() +
714  fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
715  if (std::abs(dx) > halfCarTolerance) continue;
716  nx = fPlanes[i].a;
717  ny += fPlanes[i].b;
718  nz += fPlanes[i].c;
719  ++nsurf;
720  break;
721  }
722  break;
723  }
724  case 2: // YZ section is a rectangle and
725  { // XZ section is an isosceles trapezoid
726  G4double dy = std::abs(p.y()) + fPlanes[1].d;
727  if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
728  G4double dx = fPlanes[3].a*std::abs(p.x()) +
729  fPlanes[3].c*p.z() + fPlanes[3].d;
730  if (std::abs(dx) <= halfCarTolerance)
731  {
732  nx = (p.x() < 0) ? -fPlanes[3].a : fPlanes[3].a;
733  nz += fPlanes[3].c;
734  ++nsurf;
735  }
736  break;
737  }
738  case 3: // YZ section is a rectangle and
739  { // XY section is an isosceles trapezoid
740  G4double dy = std::abs(p.y()) + fPlanes[1].d;
741  if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
742  G4double dx = fPlanes[3].a*std::abs(p.x()) +
743  fPlanes[3].b*p.y() + fPlanes[3].d;
744  if (std::abs(dx) <= halfCarTolerance)
745  {
746  nx = (p.x() < 0) ? -fPlanes[3].a : fPlanes[3].a;
747  ny += fPlanes[3].b;
748  ++nsurf;
749  }
750  break;
751  }
752  }
753 
754  // Return normal
755  //
756  if (nsurf == 1) return G4ThreeVector(nx,ny,nz);
757  else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
758  else
759  {
760  // Point is not on the surface
761  //
762 #ifdef G4CSGDEBUG
763  std::ostringstream message;
764  G4int oldprc = message.precision(16);
765  message << "Point p is not on surface (!?) of solid: "
766  << GetName() << G4endl;
767  message << "Position:\n";
768  message << " p.x() = " << p.x()/mm << " mm\n";
769  message << " p.y() = " << p.y()/mm << " mm\n";
770  message << " p.z() = " << p.z()/mm << " mm";
771  G4cout.precision(oldprc) ;
772  G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
773  JustWarning, message );
774  DumpInfo();
775 #endif
776  return ApproxSurfaceNormal(p);
777  }
778 }
779 
781 //
782 // Algorithm for SurfaceNormal() following the original specification
783 // for points not on the surface
784 
786 {
787  G4double dist = -DBL_MAX;
788  G4int iside = 0;
789  for (G4int i=0; i<4; ++i)
790  {
791  G4double d = fPlanes[i].a*p.x() +
792  fPlanes[i].b*p.y() +
793  fPlanes[i].c*p.z() + fPlanes[i].d;
794  if (d > dist) { dist = d; iside = i; }
795  }
796 
797  G4double distz = std::abs(p.z()) - fDz;
798  if (dist > distz)
799  return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
800  else
801  return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
802 }
803 
805 //
806 // Calculate distance to shape from outside
807 // - return kInfinity if no intersection
808 
810  const G4ThreeVector& v ) const
811 {
812  // Z intersections
813  //
814  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
815  return kInfinity;
816  G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
817  G4double dz = (invz < 0) ? fDz : -fDz;
818  G4double tzmin = (p.z() + dz)*invz;
819  G4double tzmax = (p.z() - dz)*invz;
820 
821  // Y intersections
822  //
823  G4double tymin = 0, tymax = DBL_MAX;
824  G4int i = 0;
825  for ( ; i<2; ++i)
826  {
827  G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
828  G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
829  if (dist >= -halfCarTolerance)
830  {
831  if (cosa >= 0) return kInfinity;
832  G4double tmp = -dist/cosa;
833  if (tymin < tmp) tymin = tmp;
834  }
835  else if (cosa > 0)
836  {
837  G4double tmp = -dist/cosa;
838  if (tymax > tmp) tymax = tmp;
839  }
840  }
841 
842  // Z intersections
843  //
844  G4double txmin = 0, txmax = DBL_MAX;
845  for ( ; i<4; ++i)
846  {
847  G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
848  G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() +
849  fPlanes[i].d;
850  if (dist >= -halfCarTolerance)
851  {
852  if (cosa >= 0) return kInfinity;
853  G4double tmp = -dist/cosa;
854  if (txmin < tmp) txmin = tmp;
855  }
856  else if (cosa > 0)
857  {
858  G4double tmp = -dist/cosa;
859  if (txmax > tmp) txmax = tmp;
860  }
861  }
862 
863  // Find distance
864  //
865  G4double tmin = std::max(std::max(txmin,tymin),tzmin);
866  G4double tmax = std::min(std::min(txmax,tymax),tzmax);
867 
868  if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
869  return (tmin < halfCarTolerance ) ? 0. : tmin;
870 }
871 
873 //
874 // Calculate exact shortest distance to any boundary from outside
875 // This is the best fast estimation of the shortest distance to trap
876 // - return 0 if point is inside
877 
879 {
880  switch (fTrapType)
881  {
882  case 0: // General case
883  {
884  G4double dz = std::abs(p.z())-fDz;
885  G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
886  G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
887  G4double dy = std::max(dz,std::max(dy1,dy2));
888 
889  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
890  + fPlanes[2].c*p.z()+fPlanes[2].d;
891  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
892  + fPlanes[3].c*p.z()+fPlanes[3].d;
893  G4double dist = std::max(dy,std::max(dx1,dx2));
894  return (dist > 0) ? dist : 0.;
895  }
896  case 1: // YZ section is a rectangle
897  {
898  G4double dz = std::abs(p.z())-fDz;
899  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
900  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
901  + fPlanes[2].c*p.z()+fPlanes[2].d;
902  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
903  + fPlanes[3].c*p.z()+fPlanes[3].d;
904  G4double dist = std::max(dy,std::max(dx1,dx2));
905  return (dist > 0) ? dist : 0.;
906  }
907  case 2: // YZ section is a rectangle and
908  { // XZ section is an isosceles trapezoid
909  G4double dz = std::abs(p.z())-fDz;
910  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
911  G4double dx = fPlanes[3].a*std::abs(p.x())
912  + fPlanes[3].c*p.z()+fPlanes[3].d;
913  G4double dist = std::max(dy,dx);
914  return (dist > 0) ? dist : 0.;
915  }
916  case 3: // YZ section is a rectangle and
917  { // XY section is an isosceles trapezoid
918  G4double dz = std::abs(p.z())-fDz;
919  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
920  G4double dx = fPlanes[3].a*std::abs(p.x())
921  + fPlanes[3].b*p.y()+fPlanes[3].d;
922  G4double dist = std::max(dy,dx);
923  return (dist > 0) ? dist : 0.;
924  }
925  }
926  return 0.;
927 }
928 
930 //
931 // Calculate distance to surface of shape from inside and
932 // find normal at exit point, if required
933 // - when leaving the surface, return 0
934 
936  const G4bool calcNorm,
937  G4bool* validNorm, G4ThreeVector* n) const
938 {
939  // Z intersections
940  //
941  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
942  {
943  if (calcNorm)
944  {
945  *validNorm = true;
946  n->set(0, 0, (p.z() < 0) ? -1 : 1);
947  }
948  return 0;
949  }
950  G4double vz = v.z();
951  G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
952  G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
953 
954  // Y intersections
955  //
956  G4int i = 0;
957  for ( ; i<2; ++i)
958  {
959  G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
960  if (cosa > 0)
961  {
962  G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
963  if (dist >= -halfCarTolerance)
964  {
965  if (calcNorm)
966  {
967  *validNorm = true;
968  n->set(0, fPlanes[i].b, fPlanes[i].c);
969  }
970  return 0;
971  }
972  G4double tmp = -dist/cosa;
973  if (tmax > tmp) { tmax = tmp; iside = i; }
974  }
975  }
976 
977  // X intersections
978  //
979  for ( ; i<4; ++i)
980  {
981  G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
982  if (cosa > 0)
983  {
984  G4double dist = fPlanes[i].a*p.x() +
985  fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
986  if (dist >= -halfCarTolerance)
987  {
988  if (calcNorm)
989  {
990  *validNorm = true;
991  n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
992  }
993  return 0;
994  }
995  G4double tmp = -dist/cosa;
996  if (tmax > tmp) { tmax = tmp; iside = i; }
997  }
998  }
999 
1000  // Set normal, if required, and return distance
1001  //
1002  if (calcNorm)
1003  {
1004  *validNorm = true;
1005  if (iside < 0)
1006  n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
1007  else
1008  n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
1009  }
1010  return tmax;
1011 }
1012 
1014 //
1015 // Calculate exact shortest distance to any boundary from inside
1016 // - Returns 0 is ThreeVector outside
1017 
1019 {
1020 #ifdef G4CSGDEBUG
1021  if( Inside(p) == kOutside )
1022  {
1023  std::ostringstream message;
1024  G4int oldprc = message.precision(16);
1025  message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
1026  message << "Position:\n";
1027  message << " p.x() = " << p.x()/mm << " mm\n";
1028  message << " p.y() = " << p.y()/mm << " mm\n";
1029  message << " p.z() = " << p.z()/mm << " mm";
1030  G4cout.precision(oldprc) ;
1031  G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
1032  JustWarning, message );
1033  DumpInfo();
1034  }
1035 #endif
1036  switch (fTrapType)
1037  {
1038  case 0: // General case
1039  {
1040  G4double dz = std::abs(p.z())-fDz;
1041  G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1042  G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1043  G4double dy = std::max(dz,std::max(dy1,dy2));
1044 
1045  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1046  + fPlanes[2].c*p.z()+fPlanes[2].d;
1047  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1048  + fPlanes[3].c*p.z()+fPlanes[3].d;
1049  G4double dist = std::max(dy,std::max(dx1,dx2));
1050  return (dist < 0) ? -dist : 0.;
1051  }
1052  case 1: // YZ section is a rectangle
1053  {
1054  G4double dz = std::abs(p.z())-fDz;
1055  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1056  G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1057  + fPlanes[2].c*p.z()+fPlanes[2].d;
1058  G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1059  + fPlanes[3].c*p.z()+fPlanes[3].d;
1060  G4double dist = std::max(dy,std::max(dx1,dx2));
1061  return (dist < 0) ? -dist : 0.;
1062  }
1063  case 2: // YZ section is a rectangle and
1064  { // XZ section is an isosceles trapezoid
1065  G4double dz = std::abs(p.z())-fDz;
1066  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1067  G4double dx = fPlanes[3].a*std::abs(p.x())
1068  + fPlanes[3].c*p.z()+fPlanes[3].d;
1069  G4double dist = std::max(dy,dx);
1070  return (dist < 0) ? -dist : 0.;
1071  }
1072  case 3: // YZ section is a rectangle and
1073  { // XY section is an isosceles trapezoid
1074  G4double dz = std::abs(p.z())-fDz;
1075  G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1076  G4double dx = fPlanes[3].a*std::abs(p.x())
1077  + fPlanes[3].b*p.y()+fPlanes[3].d;
1078  G4double dist = std::max(dy,dx);
1079  return (dist < 0) ? -dist : 0.;
1080  }
1081  }
1082  return 0.;
1083 }
1084 
1086 //
1087 // GetEntityType
1088 
1090 {
1091  return G4String("G4Trap");
1092 }
1093 
1095 //
1096 // Make a clone of the object
1097 //
1099 {
1100  return new G4Trap(*this);
1101 }
1102 
1104 //
1105 // Stream object contents to an output stream
1106 
1107 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1108 {
1109  G4double phi = std::atan2(fTthetaSphi,fTthetaCphi);
1110  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1112  G4double alpha1 = std::atan(fTalpha1);
1113  G4double alpha2 = std::atan(fTalpha2);
1114  G4String signDegree = "\u00B0";
1115 
1116  G4int oldprc = os.precision(16);
1117  os << "-----------------------------------------------------------\n"
1118  << " *** Dump for solid: " << GetName() << " ***\n"
1119  << " ===================================================\n"
1120  << " Solid type: G4Trap\n"
1121  << " Parameters:\n"
1122  << " half length Z: " << fDz/mm << " mm\n"
1123  << " half length Y, face -Dz: " << fDy1/mm << " mm\n"
1124  << " half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n"
1125  << " half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n"
1126  << " half length Y, face +Dz: " << fDy2/mm << " mm\n"
1127  << " half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n"
1128  << " half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n"
1129  << " theta: " << theta/degree << signDegree << "\n"
1130  << " phi: " << phi/degree << signDegree << "\n"
1131  << " alpha, face -Dz: " << alpha1/degree << signDegree << "\n"
1132  << " alpha, face +Dz: " << alpha2/degree << signDegree << "\n"
1133  << "-----------------------------------------------------------\n";
1134  os.precision(oldprc);
1135 
1136  return os;
1137 }
1138 
1140 //
1141 // Compute vertices from planes
1142 
1144 {
1145  for (G4int i=0; i<8; ++i)
1146  {
1147  G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
1148  G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
1149  G4double z = (i < 4) ? -fDz : fDz;
1150  G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b;
1151  G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z
1152  + fPlanes[ix].d)/fPlanes[ix].a;
1153  pt[i].set(x,y,z);
1154  }
1155 }
1156 
1158 //
1159 // Generate random point on the surface
1160 
1162 {
1163  G4ThreeVector pt[8];
1164  G4int iface [6][4] =
1165  { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
1166  G4double sface[6];
1167 
1168  GetVertices(pt);
1169  G4double stotal = 0;
1170  for (G4int i=0; i<6; ++i)
1171  {
1172  G4double ss = G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
1173  pt[iface[i][1]],
1174  pt[iface[i][2]],
1175  pt[iface[i][3]]).mag();
1176  stotal += ss;
1177  sface[i] = stotal;
1178  }
1179 
1180  // Select face
1181  //
1182  G4double select = stotal*G4UniformRand();
1183  G4int k = 5;
1184  if (select <= sface[4]) k = 4;
1185  if (select <= sface[3]) k = 3;
1186  if (select <= sface[2]) k = 2;
1187  if (select <= sface[1]) k = 1;
1188  if (select <= sface[0]) k = 0;
1189 
1190  // Select sub-triangle
1191  //
1192  G4int i0 = iface[k][0];
1193  G4int i1 = iface[k][1];
1194  G4int i2 = iface[k][2];
1195  G4int i3 = iface[k][3];
1196  G4double s1 = G4GeomTools::TriangleAreaNormal(pt[i0],pt[i1],pt[i3]).mag();
1197  G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
1198  if ((s1+s2)*G4UniformRand() > s1) i0 = i2;
1199 
1200  // Generate point
1201  //
1202  G4double u = G4UniformRand();
1203  G4double v = G4UniformRand();
1204  if (u + v > 1.) { u = 1. - u; v = 1. - v; }
1205  return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
1206 }
1207 
1209 //
1210 // Methods for visualisation
1211 
1213 {
1214  scene.AddSolid (*this);
1215 }
1216 
1218 {
1219  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1220  G4double alpha1 = std::atan(fTalpha1);
1221  G4double alpha2 = std::atan(fTalpha2);
1222  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1224 
1225  return new G4PolyhedronTrap(fDz, theta, phi,
1226  fDy1, fDx1, fDx2, alpha1,
1227  fDy2, fDx3, fDx4, alpha2);
1228 }
1229 
1230 #endif