| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Implementation for G4UTrap wrapper class
27 //
28 // 13.09.13 G.Cosmo, CERN
29 // -------------------------------------------
30
31 #include "G4Trap.hh"
32 #include "G4UTrap.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G
35
36 #include "G4AffineTransform.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4BoundingEnvelope.hh"
39
40 using namespace CLHEP;
41
42 //////////////////////////////////////////////
43 //
44 // Constructors
45 //
46 G4UTrap::G4UTrap( const G4String& pName,
47 G4double pdz,
48 G4double pTheta, G4dou
49 G4double pdy1, G4doubl
50 G4double pAlp1,
51 G4double pdy2, G4doubl
52 G4double pAlp2 )
53 : Base_t(pName, pdz, pTheta, pPhi, pdy1, pdx
54 pAlp1, pdy2, pdx3, pdx4, pAlp2)
55 {
56 G4ThreeVector pt[8];
57 CheckParameters();
58 GetVertices(pt);
59 CheckPlanarity(pt);
60 }
61
62 G4UTrap::G4UTrap( const G4String& pName,
63 const G4ThreeVector pt[8] )
64 : Base_t(pName)
65 {
66 // Start with check of centering - the cente
67 // should cross the origin of frame
68 if ( pt[0].z() >= 0
69 || pt[0].z() != pt[1].z()
70 || pt[0].z() != pt[2].z()
71 || pt[0].z() != pt[3].z()
72
73 || pt[4].z() <= 0
74 || pt[4].z() != pt[5].z()
75 || pt[4].z() != pt[6].z()
76 || pt[4].z() != pt[7].z()
77
78 || std::abs( pt[0].z() + pt[4].z() ) >
79
80 || pt[0].y() != pt[1].y()
81 || pt[2].y() != pt[3].y()
82 || pt[4].y() != pt[5].y()
83 || pt[6].y() != pt[7].y()
84
85 || std::abs(pt[0].y()+pt[2].y()+pt[4].
86 || std::abs(pt[0].x()+pt[1].x()+pt[4].
87 pt[2].x()+pt[3].x()+pt[6].
88 {
89 std::ostringstream message;
90 message << "Invalid vertice coordinates fo
91 G4Exception("G4UTrap::G4UTrap()", "GeomSol
92 FatalException, message);
93 }
94
95 SetPlanes(pt);
96 CheckParameters();
97 CheckPlanarity(pt);
98 }
99
100 // Constructor for Right Angular Wedge from ST
101 // the UnplacedTrapezoid constructor taking fo
102 // a Trd1.
103 G4UTrap::G4UTrap( const G4String& pName,
104 G4double pZ,
105 G4double pY,
106 G4double pX, G4double
107 : Base_t(pName, 0.5*pZ, /*theta=*/0, /*phi=*
108 /*dx1=*/0.5*pX, /*dx2=*/0.5*pLTX, /
109 /*dy2=*/0.5*pY, /*dx3=*/0.5*pX, /*d
110 /*Alpha2=*/0.5*(pLTX - pX)/pY)
111 {
112 CheckParameters();
113 }
114
115 G4UTrap::G4UTrap( const G4String& pName,
116 G4double pdx1, G4doub
117 G4double pdy1, G4doub
118 G4double pdz )
119 : Base_t(pName, pdx1, pdx2, pdy1, pdy2, pdz)
120 {
121 CheckParameters();
122 }
123
124 G4UTrap::G4UTrap(const G4String& pName,
125 G4double pdx, G4double
126 G4double pAlpha, G4doub
127 : Base_t(pName, pdx, pdy, pdz, pAlpha, pThet
128 {
129 CheckParameters();
130 }
131
132 G4UTrap::G4UTrap( const G4String& pName )
133 : Base_t(pName)
134 {
135 }
136
137 //////////////////////////////////////////////
138 //
139 // Fake default constructor - sets only member
140 // for usage restri
141 //
142 G4UTrap::G4UTrap( __void__& a )
143 : Base_t(a)
144 {
145 }
146
147 //////////////////////////////////////////////
148 //
149 // Destructor
150 //
151 G4UTrap::~G4UTrap() = default;
152
153 //////////////////////////////////////////////
154 //
155 // Copy constructor
156 //
157 G4UTrap::G4UTrap(const G4UTrap& rhs)
158 : Base_t(rhs)
159 {
160 }
161
162 //////////////////////////////////////////////
163 //
164 // Assignment operator
165 //
166 G4UTrap& G4UTrap::operator = (const G4UTrap& r
167 {
168 // Check assignment to self
169 //
170 if (this == &rhs) { return *this; }
171
172 // Copy base class data
173 //
174 Base_t::operator=(rhs);
175
176 return *this;
177 }
178
179 //////////////////////////////////////////////
180 //
181 // Accessors
182 //
183 G4double G4UTrap::GetZHalfLength() const
184 {
185 return GetDz();
186 }
187 G4double G4UTrap::GetYHalfLength1() const
188 {
189 return GetDy1();
190 }
191 G4double G4UTrap::GetXHalfLength1() const
192 {
193 return GetDx1();
194 }
195 G4double G4UTrap::GetXHalfLength2() const
196 {
197 return GetDx2();
198 }
199 G4double G4UTrap::GetTanAlpha1() const
200 {
201 return Base_t::GetTanAlpha1();
202 }
203 G4double G4UTrap::GetYHalfLength2() const
204 {
205 return GetDy2();
206 }
207 G4double G4UTrap::GetXHalfLength3() const
208 {
209 return GetDx3();
210 }
211 G4double G4UTrap::GetXHalfLength4() const
212 {
213 return GetDx4();
214 }
215 G4double G4UTrap::GetTanAlpha2() const
216 {
217 return Base_t::GetTanAlpha2();
218 }
219 G4double G4UTrap::GetPhi() const
220 {
221 return Base_t::GetPhi();
222 }
223 G4double G4UTrap::GetTheta() const
224 {
225 return Base_t::GetTheta();
226 }
227 G4double G4UTrap::GetAlpha1() const
228 {
229 return Base_t::GetAlpha1();
230 }
231 G4double G4UTrap::GetAlpha2() const
232 {
233 return Base_t::GetAlpha2();
234 }
235 TrapSidePlane G4UTrap::GetSidePlane(G4int n) c
236 {
237 TrapSidePlane plane;
238 plane.a = GetStruct().GetPlane(n).fA;
239 plane.b = GetStruct().GetPlane(n).fB;
240 plane.c = GetStruct().GetPlane(n).fC;
241 plane.d = GetStruct().GetPlane(n).fD;
242 return plane;
243 }
244 G4ThreeVector G4UTrap::GetSymAxis() const
245 {
246 G4double tanThetaSphi = GetTanThetaSinPhi();
247 G4double tanThetaCphi = GetTanThetaCosPhi();
248 G4double tan2Theta = tanThetaSphi*tanThetaSp
249 G4double cosTheta = 1.0 / std::sqrt(1 + tan2
250 return {tanThetaCphi*cosTheta, tanThetaSphi*
251 }
252
253 //////////////////////////////////////////////
254 //
255 // Modifier
256 //
257 void G4UTrap::SetAllParameters(G4double pDz, G
258 G4double pDy1,
259 G4double pAlp1,
260 G4double pDy2,
261 G4double pAlp2)
262 {
263 SetDz(pDz);
264 SetDy1(pDy1);
265 SetDy2(pDy2);
266 SetDx1(pDx1);
267 SetDx2(pDx2);
268 SetDx3(pDx3);
269 SetDx4(pDx4);
270 SetTanAlpha1(std::tan(pAlp1));
271 SetTanAlpha2(std::tan(pAlp2));
272 // last two will also reset cached variables
273 SetTheta(pTheta);
274 SetPhi(pPhi);
275 fRebuildPolyhedron = true;
276
277 G4ThreeVector pt[8];
278 CheckParameters();
279 GetVertices(pt);
280 CheckPlanarity(pt);
281 }
282
283 //////////////////////////////////////////////
284 //
285 // Set parameters using eight vertices
286 //
287 void G4UTrap::SetPlanes(const G4ThreeVector pt
288 {
289 U3Vector upt[8];
290 for (unsigned int i=0; i<8; ++i)
291 {
292 upt[i] = U3Vector(pt[i].x(), pt[i].y(), pt
293 }
294 fromCornersToParameters(upt);
295 fRebuildPolyhedron = true;
296 }
297
298 //////////////////////////////////////////////
299 //
300 // Check dimensions
301 //
302 void G4UTrap::CheckParameters() const
303 {
304 G4double fDz = GetZHalfLength();
305 G4double fDy1 = GetYHalfLength1();
306 G4double fDx1 = GetXHalfLength1();
307 G4double fDx2 = GetXHalfLength2();
308 G4double fDy2 = GetYHalfLength2();
309 G4double fDx3 = GetXHalfLength3();
310 G4double fDx4 = GetXHalfLength4();
311
312 if (fDz<=0 ||
313 fDy1<=0 || fDx1<=0 || fDx2<=0 ||
314 fDy2<=0 || fDx3<=0 || fDx4<=0)
315 {
316 std::ostringstream message;
317 message << "Invalid Length Parameters for
318 << "\n X - " <<fDx1<<", "<<fDx2<<
319 << "\n Y - " <<fDy1<<", "<<fDy2
320 << "\n Z - " <<fDz;
321 G4Exception("G4UTrap::CheckParameters()",
322 FatalException, message);
323 }
324 }
325
326 //////////////////////////////////////////////
327 //
328 // Compute coordinates of vertices
329 //
330 void G4UTrap::GetVertices(G4ThreeVector pt[8])
331 {
332 G4double fDz = GetZHalfLength();
333 G4double fDy1 = GetYHalfLength1();
334 G4double fDx1 = GetXHalfLength1();
335 G4double fDx2 = GetXHalfLength2();
336 G4double fDy2 = GetYHalfLength2();
337 G4double fDx3 = GetXHalfLength3();
338 G4double fDx4 = GetXHalfLength4();
339 G4double fTalpha1 = GetTanAlpha1();
340 G4double fTalpha2 = GetTanAlpha2();
341
342 G4double DzTthetaCphi = fDz*GetTanThetaCosPh
343 G4double DzTthetaSphi = fDz*GetTanThetaSinPh
344 G4double Dy1Talpha1 = fDy1*fTalpha1;
345 G4double Dy2Talpha2 = fDy2*fTalpha2;
346
347 pt[0].set(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzT
348 pt[1].set(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzT
349 pt[2].set(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzT
350 pt[3].set(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzT
351 pt[4].set( DzTthetaCphi-Dy2Talpha2-fDx3, DzT
352 pt[5].set( DzTthetaCphi-Dy2Talpha2+fDx3, DzT
353 pt[6].set( DzTthetaCphi+Dy2Talpha2-fDx4, DzT
354 pt[7].set( DzTthetaCphi+Dy2Talpha2+fDx4, DzT
355 }
356
357 //////////////////////////////////////////////
358 //
359 // Check planarity of lateral planes
360 //
361 void G4UTrap::CheckPlanarity(const G4ThreeVect
362 {
363 constexpr G4int iface[4][4] = { {0,4,5,1}, {
364 const static G4String side[4] = { "~-Y", "~+
365
366 for (G4int i=0; i<4; ++i)
367 {
368 TrapSidePlane plane = GetSidePlane(i);
369 G4double dmax = 0;
370 for (G4int k=0; k<4; ++k)
371 {
372 const G4ThreeVector p = pt[iface[i][k]];
373 G4double dist = plane.a*p.x() + plane.b*
374 if (std::abs(dist) > std::abs(dmax)) dma
375 }
376 if (std::abs(dmax) > 1000 * kCarTolerance)
377 {
378 std::ostringstream message;
379 message << "Side face " << side[i] << "
380 << GetName() << "\nDiscrepancy:
381 StreamInfo(message);
382 G4Exception("G4UTrap::CheckPlanarity()",
383 FatalException, message);
384 }
385 }
386 }
387
388 //////////////////////////////////////////////
389 //
390 // Dispatch to parameterisation for replicatio
391 // computation & modification.
392 //
393 void G4UTrap::ComputeDimensions( G4VPVPar
394 const G4int n,
395 const G4VPhysi
396 {
397 p->ComputeDimensions(*(G4Trap*)this,n,pRep);
398 }
399
400 //////////////////////////////////////////////
401 //
402 // Make a clone of the object
403 //
404 G4VSolid* G4UTrap::Clone() const
405 {
406 return new G4UTrap(*this);
407 }
408
409 //////////////////////////////////////////////
410 //
411 // Get bounding box
412
413 void G4UTrap::BoundingLimits(G4ThreeVector& pM
414 {
415 static G4bool checkBBox = true;
416
417 TrapSidePlane planes[4];
418 for (G4int i=0; i<4; ++i) { planes[i] = GetS
419
420 G4double xmin = kInfinity, xmax = -kInfinity
421 G4double ymin = kInfinity, ymax = -kInfinity
422 G4double dz = GetZHalfLength();
423 for (G4int i=0; i<8; ++i)
424 {
425 G4int iy = (i==0 || i==1 || i==4 || i==5)
426 G4int ix = (i==0 || i==2 || i==4 || i==6)
427 G4double z = (i < 4) ? -dz : dz;
428 G4double y = -(planes[iy].c*z + planes[iy]
429 G4double x = -(planes[ix].b*y + planes[ix]
430 if (x < xmin) xmin = x;
431 if (x > xmax) xmax = x;
432 if (y < ymin) ymin = y;
433 if (y > ymax) ymax = y;
434 }
435
436 pMin.set(xmin,ymin,-dz);
437 pMax.set(xmax,ymax, dz);
438
439 // Check correctness of the bounding box
440 //
441 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
442 {
443 std::ostringstream message;
444 message << "Bad bounding box (min >= max)
445 << GetName() << " !"
446 << "\npMin = " << pMin
447 << "\npMax = " << pMax;
448 G4Exception("G4UTrap::BoundingLimits()", "
449 JustWarning, message);
450 StreamInfo(G4cout);
451 }
452
453 // Check consistency of bounding boxes
454 //
455 if (checkBBox)
456 {
457 G4double tolerance = kCarTolerance;
458 U3Vector vmin, vmax;
459 Extent(vmin,vmax);
460 if (std::abs(pMin.x()-vmin.x()) > toleranc
461 std::abs(pMin.y()-vmin.y()) > toleranc
462 std::abs(pMin.z()-vmin.z()) > toleranc
463 std::abs(pMax.x()-vmax.x()) > toleranc
464 std::abs(pMax.y()-vmax.y()) > toleranc
465 std::abs(pMax.z()-vmax.z()) > toleranc
466 {
467 std::ostringstream message;
468 message << "Inconsistency in bounding bo
469 << GetName() << " !"
470 << "\nBBox min: wrapper = " << p
471 << "\nBBox max: wrapper = " << p
472 G4Exception("G4UTrap::BoundingLimits()",
473 JustWarning, message);
474 checkBBox = false;
475 }
476 }
477 }
478
479 //////////////////////////////////////////////
480 //
481 // Calculate extent under transform and specif
482
483 G4bool
484 G4UTrap::CalculateExtent(const EAxis pAxis,
485 const G4VoxelLimits&
486 const G4AffineTransfo
487 G4double& pMin,
488 {
489 G4ThreeVector bmin, bmax;
490 G4bool exist;
491
492 // Check bounding box (bbox)
493 //
494 BoundingLimits(bmin,bmax);
495 G4BoundingEnvelope bbox(bmin,bmax);
496 #ifdef G4BBOX_EXTENT
497 if (true) return bbox.CalculateExtent(pAxis,
498 #endif
499 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
500 {
501 return exist = pMin < pMax;
502 }
503
504 // Set bounding envelope (benv) and calculat
505 //
506 TrapSidePlane planes[4];
507 for (G4int i=0; i<4; ++i) { planes[i] = GetS
508
509 G4ThreeVector pt[8];
510 G4double dz = GetZHalfLength();
511 for (G4int i=0; i<8; ++i)
512 {
513 G4int iy = (i==0 || i==1 || i==4 || i==5)
514 G4int ix = (i==0 || i==2 || i==4 || i==6)
515 G4double z = (i < 4) ? -dz : dz;
516 G4double y = -(planes[iy].c*z + planes[iy]
517 G4double x = -(planes[ix].b*y + planes[ix]
518 pt[i].set(x,y,z);
519 }
520
521 G4ThreeVectorList baseA(4), baseB(4);
522 baseA[0] = pt[0];
523 baseA[1] = pt[1];
524 baseA[2] = pt[3];
525 baseA[3] = pt[2];
526
527 baseB[0] = pt[4];
528 baseB[1] = pt[5];
529 baseB[2] = pt[7];
530 baseB[3] = pt[6];
531
532 std::vector<const G4ThreeVectorList *> polyg
533 polygons[0] = &baseA;
534 polygons[1] = &baseB;
535
536 G4BoundingEnvelope benv(bmin,bmax,polygons);
537 exist = benv.CalculateExtent(pAxis,pVoxelLim
538 return exist;
539 }
540
541 //////////////////////////////////////////////
542 //
543 // Create polyhedron for visualization
544 //
545 G4Polyhedron* G4UTrap::CreatePolyhedron() cons
546 {
547 return new G4PolyhedronTrap(GetZHalfLength()
548 GetYHalfLength1(
549 GetXHalfLength1(
550 GetYHalfLength2(
551 GetXHalfLength3(
552 }
553
554 #endif // G4GEOM_USE_USOLIDS
555