| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation for G4UPara wrapper class
27 //
28 // 13.09.13 G.Cosmo, CERN/PH
29 // --------------------------------------------------------------------
30
31 #include "G4Para.hh"
32 #include "G4UPara.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35
36 #include "G4AffineTransform.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4BoundingEnvelope.hh"
39
40 using namespace CLHEP;
41
42 //////////////////////////////////////////////////////////////////////////
43 //
44 // Constructor - set & check half widths
45
46 G4UPara::G4UPara(const G4String& pName,
47 G4double pDx, G4double pDy, G4double pDz,
48 G4double pAlpha, G4double pTheta, G4double pPhi)
49 : Base_t(pName, pDx, pDy, pDz, pAlpha, pTheta, pPhi)
50 {
51 fTalpha = std::tan(pAlpha);
52 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
53 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
54 CheckParameters();
55 MakePlanes();
56 }
57
58 //////////////////////////////////////////////////////////////////////////
59 //
60 // Constructor - design of trapezoid based on 8 vertices
61
62 G4UPara::G4UPara( const G4String& pName,
63 const G4ThreeVector pt[8] )
64 : Base_t(pName)
65 {
66 // Find dimensions and trigonometric values
67 //
68 G4double fDx = (pt[3].x() - pt[2].x())*0.5;
69 G4double fDy = (pt[2].y() - pt[1].y())*0.5;
70 G4double fDz = pt[7].z();
71 SetDimensions(fDx, fDy, fDz);
72 CheckParameters(); // check dimensions
73
74 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() - pt[0].x())*0.25/fDy;
75 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx)/fDz;
76 fTthetaSphi = (pt[4].y() + fDy)/fDz;
77 SetAlpha(std::atan(fTalpha));
78 SetTheta(std::atan(std::sqrt(fTthetaSphi*fTthetaSphi
79 + fTthetaCphi*fTthetaCphi)));
80 SetPhi (std::atan2(fTthetaSphi, fTthetaCphi));
81 MakePlanes();
82
83 // Recompute vertices
84 //
85 G4ThreeVector v[8];
86 G4double DyTalpha = fDy*fTalpha;
87 G4double DzTthetaSphi = fDz*fTthetaSphi;
88 G4double DzTthetaCphi = fDz*fTthetaCphi;
89 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
90 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
91 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
92 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
93 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz);
94 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz);
95 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz);
96 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz);
97
98 // Compare with original vertices
99 //
100 for (G4int i=0; i<8; ++i)
101 {
102 G4double delx = std::abs(pt[i].x() - v[i].x());
103 G4double dely = std::abs(pt[i].y() - v[i].y());
104 G4double delz = std::abs(pt[i].z() - v[i].z());
105 G4double discrepancy = std::max(std::max(delx,dely),delz);
106 if (discrepancy > 0.1*kCarTolerance)
107 {
108 std::ostringstream message;
109 G4long oldprc = message.precision(16);
110 message << "Invalid vertice coordinates for Solid: " << GetName()
111 << "\nVertix #" << i << ", discrepancy = " << discrepancy
112 << "\n original : " << pt[i]
113 << "\n recomputed : " << v[i];
114 G4cout.precision(oldprc);
115 G4Exception("G4UPara::G4UPara()", "GeomSolids0002",
116 FatalException, message);
117
118 }
119 }
120 }
121
122 //////////////////////////////////////////////////////////////////////////
123 //
124 // Fake default constructor - sets only member data and allocates memory
125 // for usage restricted to object persistency
126
127 G4UPara::G4UPara( __void__& a )
128 : Base_t(a)
129 {
130 SetAllParameters(1., 1., 1., 0., 0., 0.);
131 fRebuildPolyhedron = false;
132 }
133
134 //////////////////////////////////////////////////////////////////////////
135 //
136 // Destructor
137
138 G4UPara::~G4UPara() = default;
139
140 //////////////////////////////////////////////////////////////////////////
141 //
142 // Copy constructor
143
144 G4UPara::G4UPara(const G4UPara& rhs)
145 : Base_t(rhs), fTalpha(rhs.fTalpha),
146 fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(rhs.fTthetaSphi)
147 {
148 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
149 }
150
151 //////////////////////////////////////////////////////////////////////////
152 //
153 // Assignment operator
154
155 G4UPara& G4UPara::operator = (const G4UPara& rhs)
156 {
157 // Check assignment to self
158 //
159 if (this == &rhs) { return *this; }
160
161 // Copy base class data
162 //
163 Base_t::operator=(rhs);
164
165 // Copy data
166 //
167 fTalpha = rhs.fTalpha;
168 fTthetaCphi = rhs.fTthetaCphi;
169 fTthetaSphi = rhs.fTthetaSphi;
170 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
171
172 return *this;
173 }
174
175 //////////////////////////////////////////////////////////////////////////
176 //
177 // Accessors & modifiers
178
179 G4double G4UPara::GetZHalfLength() const
180 {
181 return GetZ();
182 }
183 G4double G4UPara::GetYHalfLength() const
184 {
185 return GetY();
186 }
187 G4double G4UPara::GetXHalfLength() const
188 {
189 return GetX();
190 }
191 G4ThreeVector G4UPara::GetSymAxis() const
192 {
193 return G4ThreeVector(fTthetaCphi,fTthetaSphi,1.).unit();
194 }
195 G4double G4UPara::GetTanAlpha() const
196 {
197 return fTalpha;
198 }
199
200 G4double G4UPara::GetPhi() const
201 {
202 return std::atan2(fTthetaSphi,fTthetaCphi);
203 }
204
205 G4double G4UPara::GetTheta() const
206 {
207 return std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
208 +fTthetaSphi*fTthetaSphi));
209 }
210
211 G4double G4UPara::GetAlpha() const
212 {
213 return std::atan(fTalpha);
214 }
215
216 void G4UPara::SetXHalfLength(G4double val)
217 {
218 SetDimensions(val, GetY(), GetZ());
219 fRebuildPolyhedron = true;
220
221 CheckParameters();
222 MakePlanes();
223 }
224 void G4UPara::SetYHalfLength(G4double val)
225 {
226 SetDimensions(GetX(), val, GetZ());
227 fRebuildPolyhedron = true;
228
229 CheckParameters();
230 MakePlanes();
231 }
232 void G4UPara::SetZHalfLength(G4double val)
233 {
234 SetDimensions(GetX(), GetY(), val);
235 fRebuildPolyhedron = true;
236
237 CheckParameters();
238 MakePlanes();
239 }
240 void G4UPara::SetAlpha(G4double alpha)
241 {
242 Base_t::SetAlpha(alpha);
243 fTalpha = std::tan(alpha);
244 fRebuildPolyhedron = true;
245
246 MakePlanes();
247 }
248 void G4UPara::SetTanAlpha(G4double val)
249 {
250 fTalpha = val;
251 fRebuildPolyhedron = true;
252
253 MakePlanes();
254 }
255 void G4UPara::SetThetaAndPhi(double pTheta, double pPhi)
256 {
257 Base_t::SetThetaAndPhi(pTheta, pPhi);
258 G4double tanTheta = std::tan(pTheta);
259 fTthetaCphi = tanTheta*std::cos(pPhi);
260 fTthetaSphi = tanTheta*std::sin(pPhi);
261 fRebuildPolyhedron = true;
262
263 MakePlanes();
264 }
265
266 //////////////////////////////////////////////////////////////////////////
267 //
268 // Set all parameters, as for constructor - set and check half-widths
269
270 void G4UPara::SetAllParameters(G4double pDx, G4double pDy, G4double pDz,
271 G4double pAlpha, G4double pTheta, G4double pPhi)
272 {
273 // Reset data of the base class
274 fRebuildPolyhedron = true;
275
276 // Set parameters
277 SetDimensions(pDx, pDy, pDz);
278 Base_t::SetAlpha(pAlpha);
279 Base_t::SetThetaAndPhi(pTheta, pPhi);
280 fTalpha = std::tan(pAlpha);
281 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
282 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
283
284 CheckParameters();
285 MakePlanes();
286 }
287
288 //////////////////////////////////////////////////////////////////////////
289 //
290 // Check dimensions
291
292 void G4UPara::CheckParameters()
293 {
294 if (GetX() < 2*kCarTolerance ||
295 GetY() < 2*kCarTolerance ||
296 GetZ() < 2*kCarTolerance)
297 {
298 std::ostringstream message;
299 message << "Invalid (too small or negative) dimensions for Solid: "
300 << GetName()
301 << "\n X - " << GetX()
302 << "\n Y - " << GetY()
303 << "\n Z - " << GetZ();
304 G4Exception("G4UPara::CheckParameters()", "GeomSolids0002",
305 FatalException, message);
306 }
307 }
308
309 //////////////////////////////////////////////////////////////////////////
310 //
311 // Set side planes
312
313 void G4UPara::MakePlanes()
314 {
315 G4ThreeVector vx(1, 0, 0);
316 G4ThreeVector vy(fTalpha, 1, 0);
317 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1);
318
319 // Set -Y & +Y planes
320 //
321 G4ThreeVector ynorm = (vx.cross(vz)).unit();
322
323 fPlanes[0].a = 0.;
324 fPlanes[0].b = ynorm.y();
325 fPlanes[0].c = ynorm.z();
326 fPlanes[0].d = fPlanes[0].b*GetY(); // point (0,fDy,0) is on plane
327
328 fPlanes[1].a = 0.;
329 fPlanes[1].b = -fPlanes[0].b;
330 fPlanes[1].c = -fPlanes[0].c;
331 fPlanes[1].d = fPlanes[0].d;
332
333 // Set -X & +X planes
334 //
335 G4ThreeVector xnorm = (vz.cross(vy)).unit();
336
337 fPlanes[2].a = xnorm.x();
338 fPlanes[2].b = xnorm.y();
339 fPlanes[2].c = xnorm.z();
340 fPlanes[2].d = fPlanes[2].a*GetZ(); // point (fDx,0,0) is on plane
341
342 fPlanes[3].a = -fPlanes[2].a;
343 fPlanes[3].b = -fPlanes[2].b;
344 fPlanes[3].c = -fPlanes[2].c;
345 fPlanes[3].d = fPlanes[2].d;
346 }
347
348 //////////////////////////////////////////////////////////////////////////
349 //
350 // Dispatch to parameterisation for replication mechanism dimension
351 // computation & modification
352
353 void G4UPara::ComputeDimensions( G4VPVParameterisation* p,
354 const G4int n,
355 const G4VPhysicalVolume* pRep )
356 {
357 p->ComputeDimensions(*(G4Para*)this,n,pRep);
358 }
359
360 //////////////////////////////////////////////////////////////////////////
361 //
362 // Get bounding box
363
364 void G4UPara::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
365 {
366 G4double dz = GetZHalfLength();
367 G4double dx = GetXHalfLength();
368 G4double dy = GetYHalfLength();
369
370 G4double x0 = dz*fTthetaCphi;
371 G4double x1 = dy*GetTanAlpha();
372 G4double xmin =
373 std::min(
374 std::min(
375 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx);
376 G4double xmax =
377 std::max(
378 std::max(
379 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx);
380
381 G4double y0 = dz*fTthetaSphi;
382 G4double ymin = std::min(-y0-dy,y0-dy);
383 G4double ymax = std::max(-y0+dy,y0+dy);
384
385 pMin.set(xmin,ymin,-dz);
386 pMax.set(xmax,ymax, dz);
387
388 // Check correctness of the bounding box
389 //
390 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
391 {
392 std::ostringstream message;
393 message << "Bad bounding box (min >= max) for solid: "
394 << GetName() << " !"
395 << "\npMin = " << pMin
396 << "\npMax = " << pMax;
397 G4Exception("G4UPara::BoundingLimits()", "GeomMgt0001",
398 JustWarning, message);
399 StreamInfo(G4cout);
400 }
401 }
402
403 //////////////////////////////////////////////////////////////////////////
404 //
405 // Calculate extent under transform and specified limit
406
407 G4bool G4UPara::CalculateExtent( const EAxis pAxis,
408 const G4VoxelLimits& pVoxelLimit,
409 const G4AffineTransform& pTransform,
410 G4double& pMin, G4double& pMax ) const
411 {
412 G4ThreeVector bmin, bmax;
413 G4bool exist;
414
415 // Check bounding box (bbox)
416 //
417 BoundingLimits(bmin,bmax);
418 G4BoundingEnvelope bbox(bmin,bmax);
419 #ifdef G4BBOX_EXTENT
420 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
421 #endif
422 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
423 {
424 return exist = pMin < pMax;
425 }
426
427 // Set bounding envelope (benv) and calculate extent
428 //
429 G4double dz = GetZHalfLength();
430 G4double dx = GetXHalfLength();
431 G4double dy = GetYHalfLength();
432
433 G4double x0 = dz*fTthetaCphi;
434 G4double x1 = dy*GetTanAlpha();
435 G4double y0 = dz*fTthetaSphi;
436
437 G4ThreeVectorList baseA(4), baseB(4);
438 baseA[0].set(-x0-x1-dx,-y0-dy,-dz);
439 baseA[1].set(-x0-x1+dx,-y0-dy,-dz);
440 baseA[2].set(-x0+x1+dx,-y0+dy,-dz);
441 baseA[3].set(-x0+x1-dx,-y0+dy,-dz);
442
443 baseB[0].set(+x0-x1-dx, y0-dy, dz);
444 baseB[1].set(+x0-x1+dx, y0-dy, dz);
445 baseB[2].set(+x0+x1+dx, y0+dy, dz);
446 baseB[3].set(+x0+x1-dx, y0+dy, dz);
447
448 std::vector<const G4ThreeVectorList *> polygons(2);
449 polygons[0] = &baseA;
450 polygons[1] = &baseB;
451
452 G4BoundingEnvelope benv(bmin,bmax,polygons);
453 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
454 return exist;
455 }
456
457 //////////////////////////////////////////////////////////////////////////
458 //
459 // Make a clone of the object
460 //
461 G4VSolid* G4UPara::Clone() const
462 {
463 return new G4UPara(*this);
464 }
465
466 //////////////////////////////////////////////////////////////////////////
467 //
468 // Methods for visualisation
469
470 G4Polyhedron* G4UPara::CreatePolyhedron () const
471 {
472 return new G4PolyhedronPara(GetX(), GetY(), GetZ(),
473 GetAlpha(), GetTheta(), GetPhi());
474 }
475
476 #endif // G4GEOM_USE_USOLIDS
477