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