ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Para.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Para.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation for 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
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 
47 using namespace CLHEP;
48 
50 //
51 // Constructor - set & check half widths
52 
53 G4Para::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 
63 //
64 // Constructor - design of trapezoid based on 8 vertices
65 
66 G4Para::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  G4int 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 
122 //
123 // Fake default constructor - sets only member data and allocates memory
124 // for usage restricted to object persistency
125 
126 G4Para::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 
134 //
135 // Destructor
136 
138 {
139 }
140 
142 //
143 // Copy constructor
144 
146  : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
147  fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTalpha(rhs.fTalpha),
148  fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(rhs.fTthetaSphi)
149 {
150  for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
151 }
152 
154 //
155 // Assignment operator
156 
158 {
159  // Check assignment to self
160  //
161  if (this == &rhs) { return *this; }
162 
163  // Copy base class data
164  //
166 
167  // Copy data
168  //
170  fDx = rhs.fDx;
171  fDy = rhs.fDy;
172  fDz = rhs.fDz;
173  fTalpha = rhs.fTalpha;
174  fTthetaCphi = rhs.fTthetaCphi;
175  fTthetaSphi = rhs.fTthetaSphi;
176  for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
177 
178  return *this;
179 }
180 
182 //
183 // Set all parameters, as for constructor - set and check half-widths
184 
186  G4double pAlpha, G4double pTheta, G4double pPhi)
187 {
188  // Reset data of the base class
189  fCubicVolume = 0;
190  fSurfaceArea = 0;
191  fRebuildPolyhedron = true;
192 
193  // Set parameters
194  fDx = pDx;
195  fDy = pDy;
196  fDz = pDz;
197  fTalpha = std::tan(pAlpha);
198  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
199  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
200 
201  CheckParameters();
202  MakePlanes();
203 }
204 
206 //
207 // Check dimensions
208 
210 {
211  if (fDx < 2*kCarTolerance ||
212  fDy < 2*kCarTolerance ||
213  fDz < 2*kCarTolerance)
214  {
215  std::ostringstream message;
216  message << "Invalid (too small or negative) dimensions for Solid: "
217  << GetName()
218  << "\n X - " << fDx
219  << "\n Y - " << fDy
220  << "\n Z - " << fDz;
221  G4Exception("G4Para::CheckParameters()", "GeomSolids0002",
222  FatalException, message);
223  }
224 }
225 
227 //
228 // Set side planes
229 
231 {
232  G4ThreeVector vx(1, 0, 0);
233  G4ThreeVector vy(fTalpha, 1, 0);
235 
236  // Set -Y & +Y planes
237  //
238  G4ThreeVector ynorm = (vx.cross(vz)).unit();
239 
240  fPlanes[0].a = 0.;
241  fPlanes[0].b = ynorm.y();
242  fPlanes[0].c = ynorm.z();
243  fPlanes[0].d = fPlanes[0].b*fDy; // point (0,fDy,0) is on plane
244 
245  fPlanes[1].a = 0.;
246  fPlanes[1].b = -fPlanes[0].b;
247  fPlanes[1].c = -fPlanes[0].c;
248  fPlanes[1].d = fPlanes[0].d;
249 
250  // Set -X & +X planes
251  //
252  G4ThreeVector xnorm = (vz.cross(vy)).unit();
253 
254  fPlanes[2].a = xnorm.x();
255  fPlanes[2].b = xnorm.y();
256  fPlanes[2].c = xnorm.z();
257  fPlanes[2].d = fPlanes[2].a*fDx; // point (fDx,0,0) is on plane
258 
259  fPlanes[3].a = -fPlanes[2].a;
260  fPlanes[3].b = -fPlanes[2].b;
261  fPlanes[3].c = -fPlanes[2].c;
262  fPlanes[3].d = fPlanes[2].d;
263 }
264 
266 //
267 // Get volume
268 
270 {
271  // It is like G4Box, since para transformations keep the volume to be const
272  if (fCubicVolume == 0)
273  {
274  fCubicVolume = 8*fDx*fDy*fDz;
275  }
276  return fCubicVolume;
277 }
278 
280 //
281 // Get surface area
282 
284 {
285  if(fSurfaceArea == 0)
286  {
287  G4ThreeVector vx(fDx, 0, 0);
288  G4ThreeVector vy(fDy*fTalpha, fDy, 0);
290 
291  G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
292  G4double sxz = (vx.cross(vz)).mag();
293  G4double syz = (vy.cross(vz)).mag();
294 
295  fSurfaceArea = 8*(sxy+sxz+syz);
296  }
297  return fSurfaceArea;
298 }
299 
301 //
302 // Dispatch to parameterisation for replication mechanism dimension
303 // computation & modification
304 
306  const G4int n,
307  const G4VPhysicalVolume* pRep )
308 {
309  p->ComputeDimensions(*this,n,pRep);
310 }
311 
313 //
314 // Get bounding box
315 
317 {
321 
322  G4double x0 = dz*fTthetaCphi;
323  G4double x1 = dy*GetTanAlpha();
324  G4double xmin =
325  std::min(
326  std::min(
327  std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx);
328  G4double xmax =
329  std::max(
330  std::max(
331  std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx);
332 
333  G4double y0 = dz*fTthetaSphi;
334  G4double ymin = std::min(-y0-dy,y0-dy);
335  G4double ymax = std::max(-y0+dy,y0+dy);
336 
337  pMin.set(xmin,ymin,-dz);
338  pMax.set(xmax,ymax, dz);
339 
340  // Check correctness of the bounding box
341  //
342  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
343  {
344  std::ostringstream message;
345  message << "Bad bounding box (min >= max) for solid: "
346  << GetName() << " !"
347  << "\npMin = " << pMin
348  << "\npMax = " << pMax;
349  G4Exception("G4Para::BoundingLimits()", "GeomMgt0001",
350  JustWarning, message);
351  DumpInfo();
352  }
353 }
354 
356 //
357 // Calculate extent under transform and specified limit
358 
360  const G4VoxelLimits& pVoxelLimit,
361  const G4AffineTransform& pTransform,
362  G4double& pMin, G4double& pMax ) const
363 {
364  G4ThreeVector bmin, bmax;
365  G4bool exist;
366 
367  // Check bounding box (bbox)
368  //
369  BoundingLimits(bmin,bmax);
370  G4BoundingEnvelope bbox(bmin,bmax);
371 #ifdef G4BBOX_EXTENT
372  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
373 #endif
374  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
375  {
376  return exist = (pMin < pMax) ? true : false;
377  }
378 
379  // Set bounding envelope (benv) and calculate extent
380  //
384 
385  G4double x0 = dz*fTthetaCphi;
386  G4double x1 = dy*GetTanAlpha();
387  G4double y0 = dz*fTthetaSphi;
388 
389  G4ThreeVectorList baseA(4), baseB(4);
390  baseA[0].set(-x0-x1-dx,-y0-dy,-dz);
391  baseA[1].set(-x0-x1+dx,-y0-dy,-dz);
392  baseA[2].set(-x0+x1+dx,-y0+dy,-dz);
393  baseA[3].set(-x0+x1-dx,-y0+dy,-dz);
394 
395  baseB[0].set(+x0-x1-dx, y0-dy, dz);
396  baseB[1].set(+x0-x1+dx, y0-dy, dz);
397  baseB[2].set(+x0+x1+dx, y0+dy, dz);
398  baseB[3].set(+x0+x1-dx, y0+dy, dz);
399 
400  std::vector<const G4ThreeVectorList *> polygons(2);
401  polygons[0] = &baseA;
402  polygons[1] = &baseB;
403 
404  G4BoundingEnvelope benv(bmin,bmax,polygons);
405  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
406  return exist;
407 }
408 
410 //
411 // Determine where is point p, inside/on_surface/outside
412 //
413 
415 {
416  G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
417  G4double dx = std::abs(xx) + fPlanes[2].d;
418 
419  G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
420  G4double dy = std::abs(yy) + fPlanes[0].d;
421  G4double dxy = std::max(dx,dy);
422 
423  G4double dz = std::abs(p.z())-fDz;
424  G4double dist = std::max(dxy,dz);
425 
426  if (dist > halfCarTolerance) return kOutside;
427  return (dist > -halfCarTolerance) ? kSurface : kInside;
428 }
429 
431 //
432 // Determine side where point is, and return corresponding normal
433 
435 {
436  G4int nsurf = 0; // number of surfaces where p is placed
437 
438  // Check Z faces
439  //
440  G4double nz = 0;
441  G4double dz = std::abs(p.z()) - fDz;
442  if (std::abs(dz) <= halfCarTolerance)
443  {
444  nz = (p.z() < 0) ? -1 : 1;
445  ++nsurf;
446  }
447 
448  // Check Y faces
449  //
450  G4double ny = 0;
451  G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
452  if (std::abs(fPlanes[0].d + yy) <= halfCarTolerance)
453  {
454  ny = fPlanes[0].b;
455  nz += fPlanes[0].c;
456  ++nsurf;
457  }
458  else if (std::abs(fPlanes[1].d - yy) <= halfCarTolerance)
459  {
460  ny = fPlanes[1].b;
461  nz += fPlanes[1].c;
462  ++nsurf;
463  }
464 
465  // Check X faces
466  //
467  G4double nx = 0;
468  G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
469  if (std::abs(fPlanes[2].d + xx) <= halfCarTolerance)
470  {
471  nx = fPlanes[2].a;
472  ny += fPlanes[2].b;
473  nz += fPlanes[2].c;
474  ++nsurf;
475  }
476  else if (std::abs(fPlanes[3].d - xx) <= halfCarTolerance)
477  {
478  nx = fPlanes[3].a;
479  ny += fPlanes[3].b;
480  nz += fPlanes[3].c;
481  ++nsurf;
482  }
483 
484  // Return normal
485  //
486  if (nsurf == 1) return G4ThreeVector(nx,ny,nz);
487  else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
488  else
489  {
490  // Point is not on the surface
491  //
492 #ifdef G4CSGDEBUG
493  std::ostringstream message;
494  G4int oldprc = message.precision(16);
495  message << "Point p is not on surface (!?) of solid: "
496  << GetName() << G4endl;
497  message << "Position:\n";
498  message << " p.x() = " << p.x()/mm << " mm\n";
499  message << " p.y() = " << p.y()/mm << " mm\n";
500  message << " p.z() = " << p.z()/mm << " mm";
501  G4cout.precision(oldprc) ;
502  G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
503  JustWarning, message );
504  DumpInfo();
505 #endif
506  return ApproxSurfaceNormal(p);
507  }
508 }
509 
511 //
512 // Algorithm for SurfaceNormal() following the original specification
513 // for points not on the surface
514 
516 {
517  G4double dist = -DBL_MAX;
518  G4int iside = 0;
519  for (G4int i=0; i<4; ++i)
520  {
521  G4double d = fPlanes[i].a*p.x() +
522  fPlanes[i].b*p.y() +
523  fPlanes[i].c*p.z() + fPlanes[i].d;
524  if (d > dist) { dist = d; iside = i; }
525  }
526 
527  G4double distz = std::abs(p.z()) - fDz;
528  if (dist > distz)
529  return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
530  else
531  return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
532 }
533 
535 //
536 // Calculate distance to shape from outside
537 // - return kInfinity if no intersection
538 
540  const G4ThreeVector& v ) const
541 {
542  // Z intersections
543  //
544  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
545  return kInfinity;
546  G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
547  G4double dz = (invz < 0) ? fDz : -fDz;
548  G4double tzmin = (p.z() + dz)*invz;
549  G4double tzmax = (p.z() - dz)*invz;
550 
551  // Y intersections
552  //
553  G4double tmin0 = tzmin, tmax0 = tzmax;
554  G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
555  G4double disy = fPlanes[0].b*p.y() + fPlanes[0].c*p.z();
556  G4double dis0 = fPlanes[0].d + disy;
557  if (dis0 >= -halfCarTolerance)
558  {
559  if (cos0 >= 0) return kInfinity;
560  G4double tmp = -dis0/cos0;
561  if (tmin0 < tmp) tmin0 = tmp;
562  }
563  else if (cos0 > 0)
564  {
565  G4double tmp = -dis0/cos0;
566  if (tmax0 > tmp) tmax0 = tmp;
567  }
568 
569  G4double tmin1 = tmin0, tmax1 = tmax0;
570  G4double cos1 = -cos0;
571  G4double dis1 = fPlanes[1].d - disy;
572  if (dis1 >= -halfCarTolerance)
573  {
574  if (cos1 >= 0) return kInfinity;
575  G4double tmp = -dis1/cos1;
576  if (tmin1 < tmp) tmin1 = tmp;
577  }
578  else if (cos1 > 0)
579  {
580  G4double tmp = -dis1/cos1;
581  if (tmax1 > tmp) tmax1 = tmp;
582  }
583 
584  // X intersections
585  //
586  G4double tmin2 = tmin1, tmax2 = tmax1;
587  G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
588  G4double disx = fPlanes[2].a*p.x() + fPlanes[2].b*p.y() + fPlanes[2].c*p.z();
589  G4double dis2 = fPlanes[2].d + disx;
590  if (dis2 >= -halfCarTolerance)
591  {
592  if (cos2 >= 0) return kInfinity;
593  G4double tmp = -dis2/cos2;
594  if (tmin2 < tmp) tmin2 = tmp;
595  }
596  else if (cos2 > 0)
597  {
598  G4double tmp = -dis2/cos2;
599  if (tmax2 > tmp) tmax2 = tmp;
600  }
601 
602  G4double tmin3 = tmin2, tmax3 = tmax2;
603  G4double cos3 = -cos2;
604  G4double dis3 = fPlanes[3].d - disx;
605  if (dis3 >= -halfCarTolerance)
606  {
607  if (cos3 >= 0) return kInfinity;
608  G4double tmp = -dis3/cos3;
609  if (tmin3 < tmp) tmin3 = tmp;
610  }
611  else if (cos3 > 0)
612  {
613  G4double tmp = -dis3/cos3;
614  if (tmax3 > tmp) tmax3 = tmp;
615  }
616 
617  // Find distance
618  //
619  G4double tmin = tmin3, tmax = tmax3;
620  if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
621  return (tmin < halfCarTolerance ) ? 0. : tmin;
622 }
623 
625 //
626 // Calculate exact shortest distance to any boundary from outside
627 // - returns 0 is point inside
628 
630 {
631  G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
632  G4double dx = std::abs(xx) + fPlanes[2].d;
633 
634  G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
635  G4double dy = std::abs(yy) + fPlanes[0].d;
636  G4double dxy = std::max(dx,dy);
637 
638  G4double dz = std::abs(p.z())-fDz;
639  G4double dist = std::max(dxy,dz);
640 
641  return (dist > 0) ? dist : 0.;
642 }
643 
645 //
646 // Calculate distance to surface of shape from inside and, if required,
647 // find normal at exit point
648 // - when leaving the surface, return 0
649 
651  const G4bool calcNorm,
652  G4bool* validNorm, G4ThreeVector* n) const
653 {
654  // Z intersections
655  //
656  if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
657  {
658  if (calcNorm)
659  {
660  *validNorm = true;
661  n->set(0, 0, (p.z() < 0) ? -1 : 1);
662  }
663  return 0.;
664  }
665  G4double vz = v.z();
666  G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
667  G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
668 
669  // Y intersections
670  //
671  G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
672  if (cos0 > 0)
673  {
674  G4double dis0 = fPlanes[0].b*p.y() + fPlanes[0].c*p.z() + fPlanes[0].d;
675  if (dis0 >= -halfCarTolerance)
676  {
677  if (calcNorm)
678  {
679  *validNorm = true;
680  n->set(0, fPlanes[0].b, fPlanes[0].c);
681  }
682  return 0.;
683  }
684  G4double tmp = -dis0/cos0;
685  if (tmax > tmp) { tmax = tmp; iside = 0; }
686  }
687 
688  G4double cos1 = -cos0;
689  if (cos1 > 0)
690  {
691  G4double dis1 = fPlanes[1].b*p.y() + fPlanes[1].c*p.z() + fPlanes[1].d;
692  if (dis1 >= -halfCarTolerance)
693  {
694  if (calcNorm)
695  {
696  *validNorm = true;
697  n->set(0, fPlanes[1].b, fPlanes[1].c);
698  }
699  return 0.;
700  }
701  G4double tmp = -dis1/cos1;
702  if (tmax > tmp) { tmax = tmp; iside = 1; }
703  }
704 
705  // X intersections
706  //
707  G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
708  if (cos2 > 0)
709  {
710  G4double dis2 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
711  if (dis2 >= -halfCarTolerance)
712  {
713  if (calcNorm)
714  {
715  *validNorm = true;
716  n->set(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
717  }
718  return 0.;
719  }
720  G4double tmp = -dis2/cos2;
721  if (tmax > tmp) { tmax = tmp; iside = 2; }
722  }
723 
724  G4double cos3 = -cos2;
725  if (cos3 > 0)
726  {
727  G4double dis3 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
728  if (dis3 >= -halfCarTolerance)
729  {
730  if (calcNorm)
731  {
732  *validNorm = true;
733  n->set(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
734  }
735  return 0.;
736  }
737  G4double tmp = -dis3/cos3;
738  if (tmax > tmp) { tmax = tmp; iside = 3; }
739  }
740 
741  // Set normal, if required, and return distance
742  //
743  if (calcNorm)
744  {
745  *validNorm = true;
746  if (iside < 0)
747  n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
748  else
749  n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
750  }
751  return tmax;
752 }
753 
755 //
756 // Calculate exact shortest distance to any boundary from inside
757 // - returns 0 is point outside
758 
760 {
761 #ifdef G4CSGDEBUG
762  if( Inside(p) == kOutside )
763  {
764  std::ostringstream message;
765  G4int oldprc = message.precision(16);
766  message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
767  message << "Position:\n";
768  message << " p.x() = " << p.x()/mm << " mm\n";
769  message << " p.y() = " << p.y()/mm << " mm\n";
770  message << " p.z() = " << p.z()/mm << " mm";
771  G4cout.precision(oldprc) ;
772  G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
773  JustWarning, message );
774  DumpInfo();
775  }
776 #endif
777  G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
778  G4double dx = std::abs(xx) + fPlanes[2].d;
779 
780  G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
781  G4double dy = std::abs(yy) + fPlanes[0].d;
782  G4double dxy = std::max(dx,dy);
783 
784  G4double dz = std::abs(p.z())-fDz;
785  G4double dist = std::max(dxy,dz);
786 
787  return (dist < 0) ? -dist : 0.;
788 }
789 
791 //
792 // GetEntityType
793 
795 {
796  return G4String("G4Para");
797 }
798 
800 //
801 // Make a clone of the object
802 //
804 {
805  return new G4Para(*this);
806 }
807 
809 //
810 // Stream object contents to an output stream
811 
812 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
813 {
814  G4double alpha = std::atan(fTalpha);
815  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
817  G4double phi = std::atan2(fTthetaSphi,fTthetaCphi);
818  G4String signDegree = "\u00B0";
819 
820  G4int oldprc = os.precision(16);
821  os << "-----------------------------------------------------------\n"
822  << " *** Dump for solid - " << GetName() << " ***\n"
823  << " ===================================================\n"
824  << " Solid type: G4Para\n"
825  << " Parameters:\n"
826  << " half length X: " << fDx/mm << " mm\n"
827  << " half length Y: " << fDy/mm << " mm\n"
828  << " half length Z: " << fDz/mm << " mm\n"
829  << " alpha: " << alpha/degree << signDegree << "\n"
830  << " theta: " << theta/degree << signDegree << "\n"
831  << " phi: " << phi/degree << signDegree << "\n"
832  << "-----------------------------------------------------------\n";
833  os.precision(oldprc);
834 
835  return os;
836 }
837 
839 //
840 // Return a point randomly and uniformly selected on the solid surface
841 
843 {
844  G4double DyTalpha = fDy*fTalpha;
845  G4double DzTthetaSphi = fDz*fTthetaSphi;
846  G4double DzTthetaCphi = fDz*fTthetaCphi;
847 
848  // Set vertices
849  //
850  G4ThreeVector pt[8];
851  pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
852  pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
853  pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
854  pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
855  pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz);
856  pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz);
857  pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz);
858  pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz);
859 
860  // Set areas (-Z, -Y, +Y, -X, +X, +Z)
861  //
862  G4ThreeVector vx(fDx, 0, 0);
863  G4ThreeVector vy(DyTalpha, fDy, 0);
864  G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, fDz);
865 
866  G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
867  G4double sxz = (vx.cross(vz)).mag();
868  G4double syz = (vy.cross(vz)).mag();
869 
870  G4double sface[6] = { sxy, syz, syz, sxz, sxz, sxy };
871  for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; }
872 
873  // Select face
874  //
875  G4double select = sface[5]*G4UniformRand();
876  G4int k = 5;
877  if (select <= sface[4]) k = 4;
878  if (select <= sface[3]) k = 3;
879  if (select <= sface[2]) k = 2;
880  if (select <= sface[1]) k = 1;
881  if (select <= sface[0]) k = 0;
882 
883  // Generate point
884  //
885  G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, {0,2,4}, {1,5,3}, {4,6,5}};
888  return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1]] + v*pt[ip[k][2]];
889 }
890 
892 //
893 // Methods for visualisation
894 
896 {
897  scene.AddSolid (*this);
898 }
899 
901 {
902  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
903  G4double alpha = std::atan(fTalpha);
904  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
906 
907  return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
908 }
909 #endif