Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4Para.cc
Go to the documentation of this file.
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 G4Para class
27//
28// 21.03.95 P.Kent: Modified for `tolerant' geom
29// 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
30// 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
31// 29.05.17 E.Tcherniaev: complete revision, speed-up
32////////////////////////////////////////////////////////////////////////////
33
34#include "G4Para.hh"
35
36#if !defined(G4GEOM_USE_UPARA)
37
38#include "G4VoxelLimits.hh"
39#include "G4AffineTransform.hh"
40#include "G4BoundingEnvelope.hh"
41#include "Randomize.hh"
42
43#include "G4VPVParameterisation.hh"
44
45#include "G4VGraphicsScene.hh"
46
47using namespace CLHEP;
48
49//////////////////////////////////////////////////////////////////////////
50//
51// Constructor - set & check half widths
52
53G4Para::G4Para(const G4String& pName,
54 G4double pDx, G4double pDy, G4double pDz,
55 G4double pAlpha, G4double pTheta, G4double pPhi)
56 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
57{
58 SetAllParameters(pDx, pDy, pDz, pAlpha, pTheta, pPhi);
59 fRebuildPolyhedron = false; // default value for G4CSGSolid
60}
61
62//////////////////////////////////////////////////////////////////////////
63//
64// Constructor - design of trapezoid based on 8 vertices
65
66G4Para::G4Para( const G4String& pName,
67 const G4ThreeVector pt[8] )
68 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
69{
70 // Find dimensions and trigonometric values
71 //
72 fDx = (pt[3].x() - pt[2].x())*0.5;
73 fDy = (pt[2].y() - pt[1].y())*0.5;
74 fDz = pt[7].z();
75 CheckParameters(); // check dimensions
76
77 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() - pt[0].x())*0.25/fDy;
78 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx)/fDz;
79 fTthetaSphi = (pt[4].y() + fDy)/fDz;
80 MakePlanes();
81
82 // Recompute vertices
83 //
84 G4ThreeVector v[8];
85 G4double DyTalpha = fDy*fTalpha;
86 G4double DzTthetaSphi = fDz*fTthetaSphi;
87 G4double DzTthetaCphi = fDz*fTthetaCphi;
88 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
89 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
90 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
91 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
92 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz);
93 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz);
94 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz);
95 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz);
96
97 // Compare with original vertices
98 //
99 for (G4int i=0; i<8; ++i)
100 {
101 G4double delx = std::abs(pt[i].x() - v[i].x());
102 G4double dely = std::abs(pt[i].y() - v[i].y());
103 G4double delz = std::abs(pt[i].z() - v[i].z());
104 G4double discrepancy = std::max(std::max(delx,dely),delz);
105 if (discrepancy > 0.1*kCarTolerance)
106 {
107 std::ostringstream message;
108 G4long oldprc = message.precision(16);
109 message << "Invalid vertice coordinates for Solid: " << GetName()
110 << "\nVertix #" << i << ", discrepancy = " << discrepancy
111 << "\n original : " << pt[i]
112 << "\n recomputed : " << v[i];
113 G4cout.precision(oldprc);
114 G4Exception("G4Para::G4Para()", "GeomSolids0002",
115 FatalException, message);
116
117 }
118 }
119}
120
121//////////////////////////////////////////////////////////////////////////
122//
123// Fake default constructor - sets only member data and allocates memory
124// for usage restricted to object persistency
125
126G4Para::G4Para( __void__& a )
127 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance)
128{
129 SetAllParameters(1., 1., 1., 0., 0., 0.);
130 fRebuildPolyhedron = false; // default value for G4CSGSolid
131}
132
133//////////////////////////////////////////////////////////////////////////
134//
135// Destructor
136
137G4Para::~G4Para() = default;
138
139//////////////////////////////////////////////////////////////////////////
140//
141// Copy constructor
142
143G4Para::G4Para(const G4Para& rhs)
144 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
145 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), 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
155G4Para& G4Para::operator = (const G4Para& rhs)
156{
157 // Check assignment to self
158 //
159 if (this == &rhs) { return *this; }
160
161 // Copy base class data
162 //
163 G4CSGSolid::operator=(rhs);
164
165 // Copy data
166 //
167 halfCarTolerance = rhs.halfCarTolerance;
168 fDx = rhs.fDx;
169 fDy = rhs.fDy;
170 fDz = rhs.fDz;
171 fTalpha = rhs.fTalpha;
172 fTthetaCphi = rhs.fTthetaCphi;
173 fTthetaSphi = rhs.fTthetaSphi;
174 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
175
176 return *this;
177}
178
179//////////////////////////////////////////////////////////////////////////
180//
181// Set all parameters, as for constructor - set and check half-widths
182
183void G4Para::SetAllParameters(G4double pDx, G4double pDy, G4double pDz,
184 G4double pAlpha, G4double pTheta, G4double pPhi)
185{
186 // Reset data of the base class
187 fCubicVolume = 0;
188 fSurfaceArea = 0;
189 fRebuildPolyhedron = true;
190
191 // Set parameters
192 fDx = pDx;
193 fDy = pDy;
194 fDz = pDz;
195 fTalpha = std::tan(pAlpha);
196 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
197 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
198
199 CheckParameters();
200 MakePlanes();
201}
202
203//////////////////////////////////////////////////////////////////////////
204//
205// Check dimensions
206
207void G4Para::CheckParameters()
208{
209 if (fDx < 2*kCarTolerance ||
210 fDy < 2*kCarTolerance ||
211 fDz < 2*kCarTolerance)
212 {
213 std::ostringstream message;
214 message << "Invalid (too small or negative) dimensions for Solid: "
215 << GetName()
216 << "\n X - " << fDx
217 << "\n Y - " << fDy
218 << "\n Z - " << fDz;
219 G4Exception("G4Para::CheckParameters()", "GeomSolids0002",
220 FatalException, message);
221 }
222}
223
224//////////////////////////////////////////////////////////////////////////
225//
226// Set side planes
227
228void G4Para::MakePlanes()
229{
230 G4ThreeVector vx(1, 0, 0);
231 G4ThreeVector vy(fTalpha, 1, 0);
232 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1);
233
234 // Set -Y & +Y planes
235 //
236 G4ThreeVector ynorm = (vx.cross(vz)).unit();
237
238 fPlanes[0].a = 0.;
239 fPlanes[0].b = ynorm.y();
240 fPlanes[0].c = ynorm.z();
241 fPlanes[0].d = fPlanes[0].b*fDy; // point (0,fDy,0) is on plane
242
243 fPlanes[1].a = 0.;
244 fPlanes[1].b = -fPlanes[0].b;
245 fPlanes[1].c = -fPlanes[0].c;
246 fPlanes[1].d = fPlanes[0].d;
247
248 // Set -X & +X planes
249 //
250 G4ThreeVector xnorm = (vz.cross(vy)).unit();
251
252 fPlanes[2].a = xnorm.x();
253 fPlanes[2].b = xnorm.y();
254 fPlanes[2].c = xnorm.z();
255 fPlanes[2].d = fPlanes[2].a*fDx; // point (fDx,0,0) is on plane
256
257 fPlanes[3].a = -fPlanes[2].a;
258 fPlanes[3].b = -fPlanes[2].b;
259 fPlanes[3].c = -fPlanes[2].c;
260 fPlanes[3].d = fPlanes[2].d;
261}
262
263//////////////////////////////////////////////////////////////////////////
264//
265// Get volume
266
267G4double G4Para::GetCubicVolume()
268{
269 // It is like G4Box, since para transformations keep the volume to be const
270 if (fCubicVolume == 0)
271 {
272 fCubicVolume = 8*fDx*fDy*fDz;
273 }
274 return fCubicVolume;
275}
276
277//////////////////////////////////////////////////////////////////////////
278//
279// Get surface area
280
281G4double G4Para::GetSurfaceArea()
282{
283 if(fSurfaceArea == 0)
284 {
285 G4ThreeVector vx(fDx, 0, 0);
286 G4ThreeVector vy(fDy*fTalpha, fDy, 0);
287 G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTthetaSphi, fDz);
288
289 G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
290 G4double sxz = (vx.cross(vz)).mag();
291 G4double syz = (vy.cross(vz)).mag();
292
293 fSurfaceArea = 8*(sxy+sxz+syz);
294 }
295 return fSurfaceArea;
296}
297
298//////////////////////////////////////////////////////////////////////////
299//
300// Dispatch to parameterisation for replication mechanism dimension
301// computation & modification
302
303void G4Para::ComputeDimensions( G4VPVParameterisation* p,
304 const G4int n,
305 const G4VPhysicalVolume* pRep )
306{
307 p->ComputeDimensions(*this,n,pRep);
308}
309
310//////////////////////////////////////////////////////////////////////////
311//
312// Get bounding box
313
314void G4Para::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
315{
316 G4double dz = GetZHalfLength();
317 G4double dx = GetXHalfLength();
318 G4double dy = GetYHalfLength();
319
320 G4double x0 = dz*fTthetaCphi;
321 G4double x1 = dy*GetTanAlpha();
322 G4double xmin =
323 std::min(
324 std::min(
325 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx);
326 G4double xmax =
327 std::max(
328 std::max(
329 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx);
330
331 G4double y0 = dz*fTthetaSphi;
332 G4double ymin = std::min(-y0-dy,y0-dy);
333 G4double ymax = std::max(-y0+dy,y0+dy);
334
335 pMin.set(xmin,ymin,-dz);
336 pMax.set(xmax,ymax, dz);
337
338 // Check correctness of the bounding box
339 //
340 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
341 {
342 std::ostringstream message;
343 message << "Bad bounding box (min >= max) for solid: "
344 << GetName() << " !"
345 << "\npMin = " << pMin
346 << "\npMax = " << pMax;
347 G4Exception("G4Para::BoundingLimits()", "GeomMgt0001",
348 JustWarning, message);
349 DumpInfo();
350 }
351}
352
353//////////////////////////////////////////////////////////////////////////
354//
355// Calculate extent under transform and specified limit
356
357G4bool G4Para::CalculateExtent( const EAxis pAxis,
358 const G4VoxelLimits& pVoxelLimit,
359 const G4AffineTransform& pTransform,
360 G4double& pMin, G4double& pMax ) const
361{
362 G4ThreeVector bmin, bmax;
363 G4bool exist;
364
365 // Check bounding box (bbox)
366 //
367 BoundingLimits(bmin,bmax);
368 G4BoundingEnvelope bbox(bmin,bmax);
369#ifdef G4BBOX_EXTENT
370 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
371#endif
372 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
373 {
374 return exist = pMin < pMax;
375 }
376
377 // Set bounding envelope (benv) and calculate extent
378 //
379 G4double dz = GetZHalfLength();
380 G4double dx = GetXHalfLength();
381 G4double dy = GetYHalfLength();
382
383 G4double x0 = dz*fTthetaCphi;
384 G4double x1 = dy*GetTanAlpha();
385 G4double y0 = dz*fTthetaSphi;
386
387 G4ThreeVectorList baseA(4), baseB(4);
388 baseA[0].set(-x0-x1-dx,-y0-dy,-dz);
389 baseA[1].set(-x0-x1+dx,-y0-dy,-dz);
390 baseA[2].set(-x0+x1+dx,-y0+dy,-dz);
391 baseA[3].set(-x0+x1-dx,-y0+dy,-dz);
392
393 baseB[0].set(+x0-x1-dx, y0-dy, dz);
394 baseB[1].set(+x0-x1+dx, y0-dy, dz);
395 baseB[2].set(+x0+x1+dx, y0+dy, dz);
396 baseB[3].set(+x0+x1-dx, y0+dy, dz);
397
398 std::vector<const G4ThreeVectorList *> polygons(2);
399 polygons[0] = &baseA;
400 polygons[1] = &baseB;
401
402 G4BoundingEnvelope benv(bmin,bmax,polygons);
403 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
404 return exist;
405}
406
407//////////////////////////////////////////////////////////////////////////
408//
409// Determine where is point p, inside/on_surface/outside
410//
411
412EInside G4Para::Inside( const G4ThreeVector& p ) const
413{
414 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
415 G4double dx = std::abs(xx) + fPlanes[2].d;
416
417 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
418 G4double dy = std::abs(yy) + fPlanes[0].d;
419 G4double dxy = std::max(dx,dy);
420
421 G4double dz = std::abs(p.z())-fDz;
422 G4double dist = std::max(dxy,dz);
423
424 if (dist > halfCarTolerance) return kOutside;
425 return (dist > -halfCarTolerance) ? kSurface : kInside;
426}
427
428//////////////////////////////////////////////////////////////////////////
429//
430// Determine side where point is, and return corresponding normal
431
432G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
433{
434 G4int nsurf = 0; // number of surfaces where p is placed
435
436 // Check Z faces
437 //
438 G4double nz = 0;
439 G4double dz = std::abs(p.z()) - fDz;
440 if (std::abs(dz) <= halfCarTolerance)
441 {
442 nz = (p.z() < 0) ? -1 : 1;
443 ++nsurf;
444 }
445
446 // Check Y faces
447 //
448 G4double ny = 0;
449 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
450 if (std::abs(fPlanes[0].d + yy) <= halfCarTolerance)
451 {
452 ny = fPlanes[0].b;
453 nz += fPlanes[0].c;
454 ++nsurf;
455 }
456 else if (std::abs(fPlanes[1].d - yy) <= halfCarTolerance)
457 {
458 ny = fPlanes[1].b;
459 nz += fPlanes[1].c;
460 ++nsurf;
461 }
462
463 // Check X faces
464 //
465 G4double nx = 0;
466 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
467 if (std::abs(fPlanes[2].d + xx) <= halfCarTolerance)
468 {
469 nx = fPlanes[2].a;
470 ny += fPlanes[2].b;
471 nz += fPlanes[2].c;
472 ++nsurf;
473 }
474 else if (std::abs(fPlanes[3].d - xx) <= halfCarTolerance)
475 {
476 nx = fPlanes[3].a;
477 ny += fPlanes[3].b;
478 nz += fPlanes[3].c;
479 ++nsurf;
480 }
481
482 // Return normal
483 //
484 if (nsurf == 1) return {nx,ny,nz};
485 else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
486 else
487 {
488 // Point is not on the surface
489 //
490#ifdef G4CSGDEBUG
491 std::ostringstream message;
492 G4int oldprc = message.precision(16);
493 message << "Point p is not on surface (!?) of solid: "
494 << GetName() << G4endl;
495 message << "Position:\n";
496 message << " p.x() = " << p.x()/mm << " mm\n";
497 message << " p.y() = " << p.y()/mm << " mm\n";
498 message << " p.z() = " << p.z()/mm << " mm";
499 G4cout.precision(oldprc) ;
500 G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
501 JustWarning, message );
502 DumpInfo();
503#endif
504 return ApproxSurfaceNormal(p);
505 }
506}
507
508//////////////////////////////////////////////////////////////////////////
509//
510// Algorithm for SurfaceNormal() following the original specification
511// for points not on the surface
512
513G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
514{
515 G4double dist = -DBL_MAX;
516 G4int iside = 0;
517 for (G4int i=0; i<4; ++i)
518 {
519 G4double d = fPlanes[i].a*p.x() +
520 fPlanes[i].b*p.y() +
521 fPlanes[i].c*p.z() + fPlanes[i].d;
522 if (d > dist) { dist = d; iside = i; }
523 }
524
525 G4double distz = std::abs(p.z()) - fDz;
526 if (dist > distz)
527 return { fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c };
528 else
529 return { 0, 0, (G4double)((p.z() < 0) ? -1 : 1) };
530}
531
532//////////////////////////////////////////////////////////////////////////
533//
534// Calculate distance to shape from outside
535// - return kInfinity if no intersection
536
537G4double G4Para::DistanceToIn(const G4ThreeVector& p,
538 const G4ThreeVector& v ) const
539{
540 // Z intersections
541 //
542 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
543 return kInfinity;
544 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
545 G4double dz = (invz < 0) ? fDz : -fDz;
546 G4double tzmin = (p.z() + dz)*invz;
547 G4double tzmax = (p.z() - dz)*invz;
548
549 // Y intersections
550 //
551 G4double tmin0 = tzmin, tmax0 = tzmax;
552 G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
553 G4double disy = fPlanes[0].b*p.y() + fPlanes[0].c*p.z();
554 G4double dis0 = fPlanes[0].d + disy;
555 if (dis0 >= -halfCarTolerance)
556 {
557 if (cos0 >= 0) return kInfinity;
558 G4double tmp = -dis0/cos0;
559 if (tmin0 < tmp) tmin0 = tmp;
560 }
561 else if (cos0 > 0)
562 {
563 G4double tmp = -dis0/cos0;
564 if (tmax0 > tmp) tmax0 = tmp;
565 }
566
567 G4double tmin1 = tmin0, tmax1 = tmax0;
568 G4double cos1 = -cos0;
569 G4double dis1 = fPlanes[1].d - disy;
570 if (dis1 >= -halfCarTolerance)
571 {
572 if (cos1 >= 0) return kInfinity;
573 G4double tmp = -dis1/cos1;
574 if (tmin1 < tmp) tmin1 = tmp;
575 }
576 else if (cos1 > 0)
577 {
578 G4double tmp = -dis1/cos1;
579 if (tmax1 > tmp) tmax1 = tmp;
580 }
581
582 // X intersections
583 //
584 G4double tmin2 = tmin1, tmax2 = tmax1;
585 G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
586 G4double disx = fPlanes[2].a*p.x() + fPlanes[2].b*p.y() + fPlanes[2].c*p.z();
587 G4double dis2 = fPlanes[2].d + disx;
588 if (dis2 >= -halfCarTolerance)
589 {
590 if (cos2 >= 0) return kInfinity;
591 G4double tmp = -dis2/cos2;
592 if (tmin2 < tmp) tmin2 = tmp;
593 }
594 else if (cos2 > 0)
595 {
596 G4double tmp = -dis2/cos2;
597 if (tmax2 > tmp) tmax2 = tmp;
598 }
599
600 G4double tmin3 = tmin2, tmax3 = tmax2;
601 G4double cos3 = -cos2;
602 G4double dis3 = fPlanes[3].d - disx;
603 if (dis3 >= -halfCarTolerance)
604 {
605 if (cos3 >= 0) return kInfinity;
606 G4double tmp = -dis3/cos3;
607 if (tmin3 < tmp) tmin3 = tmp;
608 }
609 else if (cos3 > 0)
610 {
611 G4double tmp = -dis3/cos3;
612 if (tmax3 > tmp) tmax3 = tmp;
613 }
614
615 // Find distance
616 //
617 G4double tmin = tmin3, tmax = tmax3;
618 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
619 return (tmin < halfCarTolerance ) ? 0. : tmin;
620}
621
622//////////////////////////////////////////////////////////////////////////
623//
624// Calculate exact shortest distance to any boundary from outside
625// - returns 0 is point inside
626
627G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
628{
629 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
630 G4double dx = std::abs(xx) + fPlanes[2].d;
631
632 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
633 G4double dy = std::abs(yy) + fPlanes[0].d;
634 G4double dxy = std::max(dx,dy);
635
636 G4double dz = std::abs(p.z())-fDz;
637 G4double dist = std::max(dxy,dz);
638
639 return (dist > 0) ? dist : 0.;
640}
641
642//////////////////////////////////////////////////////////////////////////
643//
644// Calculate distance to surface of shape from inside and, if required,
645// find normal at exit point
646// - when leaving the surface, return 0
647
648G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
649 const G4bool calcNorm,
650 G4bool* validNorm, G4ThreeVector* n) const
651{
652 // Z intersections
653 //
654 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
655 {
656 if (calcNorm)
657 {
658 *validNorm = true;
659 n->set(0, 0, (p.z() < 0) ? -1 : 1);
660 }
661 return 0.;
662 }
663 G4double vz = v.z();
664 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
665 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
666
667 // Y intersections
668 //
669 G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
670 if (cos0 > 0)
671 {
672 G4double dis0 = fPlanes[0].b*p.y() + fPlanes[0].c*p.z() + fPlanes[0].d;
673 if (dis0 >= -halfCarTolerance)
674 {
675 if (calcNorm)
676 {
677 *validNorm = true;
678 n->set(0, fPlanes[0].b, fPlanes[0].c);
679 }
680 return 0.;
681 }
682 G4double tmp = -dis0/cos0;
683 if (tmax > tmp) { tmax = tmp; iside = 0; }
684 }
685
686 G4double cos1 = -cos0;
687 if (cos1 > 0)
688 {
689 G4double dis1 = fPlanes[1].b*p.y() + fPlanes[1].c*p.z() + fPlanes[1].d;
690 if (dis1 >= -halfCarTolerance)
691 {
692 if (calcNorm)
693 {
694 *validNorm = true;
695 n->set(0, fPlanes[1].b, fPlanes[1].c);
696 }
697 return 0.;
698 }
699 G4double tmp = -dis1/cos1;
700 if (tmax > tmp) { tmax = tmp; iside = 1; }
701 }
702
703 // X intersections
704 //
705 G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
706 if (cos2 > 0)
707 {
708 G4double dis2 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
709 if (dis2 >= -halfCarTolerance)
710 {
711 if (calcNorm)
712 {
713 *validNorm = true;
714 n->set(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
715 }
716 return 0.;
717 }
718 G4double tmp = -dis2/cos2;
719 if (tmax > tmp) { tmax = tmp; iside = 2; }
720 }
721
722 G4double cos3 = -cos2;
723 if (cos3 > 0)
724 {
725 G4double dis3 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
726 if (dis3 >= -halfCarTolerance)
727 {
728 if (calcNorm)
729 {
730 *validNorm = true;
731 n->set(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
732 }
733 return 0.;
734 }
735 G4double tmp = -dis3/cos3;
736 if (tmax > tmp) { tmax = tmp; iside = 3; }
737 }
738
739 // Set normal, if required, and return distance
740 //
741 if (calcNorm)
742 {
743 *validNorm = true;
744 if (iside < 0)
745 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
746 else
747 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
748 }
749 return tmax;
750}
751
752//////////////////////////////////////////////////////////////////////////
753//
754// Calculate exact shortest distance to any boundary from inside
755// - returns 0 is point outside
756
757G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
758{
759#ifdef G4CSGDEBUG
760 if( Inside(p) == kOutside )
761 {
762 std::ostringstream message;
763 G4int oldprc = message.precision(16);
764 message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
765 message << "Position:\n";
766 message << " p.x() = " << p.x()/mm << " mm\n";
767 message << " p.y() = " << p.y()/mm << " mm\n";
768 message << " p.z() = " << p.z()/mm << " mm";
769 G4cout.precision(oldprc) ;
770 G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
771 JustWarning, message );
772 DumpInfo();
773 }
774#endif
775 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
776 G4double dx = std::abs(xx) + fPlanes[2].d;
777
778 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
779 G4double dy = std::abs(yy) + fPlanes[0].d;
780 G4double dxy = std::max(dx,dy);
781
782 G4double dz = std::abs(p.z())-fDz;
783 G4double dist = std::max(dxy,dz);
784
785 return (dist < 0) ? -dist : 0.;
786}
787
788//////////////////////////////////////////////////////////////////////////
789//
790// GetEntityType
791
792G4GeometryType G4Para::GetEntityType() const
793{
794 return {"G4Para"};
795}
796
797//////////////////////////////////////////////////////////////////////////
798//
799// Make a clone of the object
800//
801G4VSolid* G4Para::Clone() const
802{
803 return new G4Para(*this);
804}
805
806//////////////////////////////////////////////////////////////////////////
807//
808// Stream object contents to an output stream
809
810std::ostream& G4Para::StreamInfo( std::ostream& os ) const
811{
812 G4double alpha = std::atan(fTalpha);
813 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
814 fTthetaSphi*fTthetaSphi));
815 G4double phi = std::atan2(fTthetaSphi,fTthetaCphi);
816
817 G4long oldprc = os.precision(16);
818 os << "-----------------------------------------------------------\n"
819 << " *** Dump for solid - " << GetName() << " ***\n"
820 << " ===================================================\n"
821 << " Solid type: G4Para\n"
822 << " Parameters:\n"
823 << " half length X: " << fDx/mm << " mm\n"
824 << " half length Y: " << fDy/mm << " mm\n"
825 << " half length Z: " << fDz/mm << " mm\n"
826 << " alpha: " << alpha/degree << "degrees\n"
827 << " theta: " << theta/degree << "degrees\n"
828 << " phi: " << phi/degree << "degrees\n"
829 << "-----------------------------------------------------------\n";
830 os.precision(oldprc);
831
832 return os;
833}
834
835//////////////////////////////////////////////////////////////////////////
836//
837// Return a point randomly and uniformly selected on the solid surface
838
839G4ThreeVector G4Para::GetPointOnSurface() const
840{
841 G4double DyTalpha = fDy*fTalpha;
842 G4double DzTthetaSphi = fDz*fTthetaSphi;
843 G4double DzTthetaCphi = fDz*fTthetaCphi;
844
845 // Set vertices
846 //
847 G4ThreeVector pt[8];
848 pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
849 pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
850 pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
851 pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
852 pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz);
853 pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz);
854 pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz);
855 pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz);
856
857 // Set areas (-Z, -Y, +Y, -X, +X, +Z)
858 //
859 G4ThreeVector vx(fDx, 0, 0);
860 G4ThreeVector vy(DyTalpha, fDy, 0);
861 G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, fDz);
862
863 G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
864 G4double sxz = (vx.cross(vz)).mag();
865 G4double syz = (vy.cross(vz)).mag();
866
867 G4double sface[6] = { sxy, syz, syz, sxz, sxz, sxy };
868 for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; }
869
870 // Select face
871 //
872 G4double select = sface[5]*G4UniformRand();
873 G4int k = 5;
874 if (select <= sface[4]) k = 4;
875 if (select <= sface[3]) k = 3;
876 if (select <= sface[2]) k = 2;
877 if (select <= sface[1]) k = 1;
878 if (select <= sface[0]) k = 0;
879
880 // Generate point
881 //
882 G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, {0,2,4}, {1,5,3}, {4,6,5}};
883 G4double u = G4UniformRand();
884 G4double v = G4UniformRand();
885 return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1]] + v*pt[ip[k][2]];
886}
887
888//////////////////////////////////////////////////////////////////////////
889//
890// Methods for visualisation
891
892void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
893{
894 scene.AddSolid (*this);
895}
896
897G4Polyhedron* G4Para::CreatePolyhedron () const
898{
899 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
900 G4double alpha = std::atan(fTalpha);
901 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
902 fTthetaSphi*fTthetaSphi));
903
904 return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
905}
906#endif
kCarTolerance
const G4double kCarTolerance
Definition G4BoundingEnvelope.cc:38
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition G4BoundingEnvelope.hh:54
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4long
long G4long
Definition G4Types.hh:87
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4BoundingEnvelope::BoundingBoxVsVoxelLimits
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Definition G4BoundingEnvelope.cc:182
G4BoundingEnvelope::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Definition G4BoundingEnvelope.cc:267
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition G4CSGSolid.hh:69
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition G4CSGSolid.hh:68
G4CSGSolid::fRebuildPolyhedron
G4bool fRebuildPolyhedron
Definition G4CSGSolid.hh:70
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition G4CSGSolid.cc:89
G4Para::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const override
Definition G4Para.cc:537
G4Para::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const override
Definition G4Para.cc:810
G4Para::~G4Para
~G4Para() override
G4Para::Inside
EInside Inside(const G4ThreeVector &p) const override
Definition G4Para.cc:412
G4Para::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const override
Definition G4Para.cc:357
G4Para::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
Definition G4Para.cc:648
G4Para::GetSurfaceArea
G4double GetSurfaceArea() override
Definition G4Para.cc:281
G4Para::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const override
Definition G4Para.cc:897
G4Para::GetTanAlpha
G4double GetTanAlpha() const
G4Para::GetEntityType
G4GeometryType GetEntityType() const override
Definition G4Para.cc:792
G4Para::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const override
Definition G4Para.cc:839
G4Para::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const override
Definition G4Para.cc:432
G4Para::operator=
G4Para & operator=(const G4Para &rhs)
Definition G4Para.cc:155
G4Para::G4Para
G4Para(const G4String &pName, G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition G4Para.cc:53
G4Para::Clone
G4VSolid * Clone() const override
Definition G4Para.cc:801
G4Para::SetAllParameters
void SetAllParameters(G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition G4Para.cc:183
G4Para::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
Definition G4Para.cc:314
G4Para::GetYHalfLength
G4double GetYHalfLength() const
G4Para::d
G4double d
Definition G4Para.hh:186
G4Para::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep) override
Definition G4Para.cc:303
G4Para::GetZHalfLength
G4double GetZHalfLength() const
G4Para::GetCubicVolume
G4double GetCubicVolume() override
Definition G4Para.cc:267
G4Para::a
G4double a
Definition G4Para.hh:186
G4Para::GetXHalfLength
G4double GetXHalfLength() const
G4Para::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const override
Definition G4Para.cc:892
G4Para::b
G4double b
Definition G4Para.hh:186
G4Para::c
G4double c
Definition G4Para.hh:186
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:89
G4VSolid::GetName
G4String GetName() const
G4VSolid::DumpInfo
void DumpInfo() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition G4VSolid.hh:299
EAxis
EAxis
Definition geomdefs.hh:54
EInside
EInside
Definition geomdefs.hh:67
kInside
@ kInside
Definition geomdefs.hh:70
kOutside
@ kOutside
Definition geomdefs.hh:68
kSurface
@ kSurface
Definition geomdefs.hh:69
DBL_MAX
#define DBL_MAX
Definition templates.hh:62