Geant4  10.02
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 //
27 // $Id: G4Para.cc 89721 2015-04-28 10:42:51Z gcosmo $
28 //
29 // class G4Para
30 //
31 // Implementation for G4Para class
32 //
33 // History:
34 //
35 // 23.10.05 V.Grichine: bug fixed in DistanceToOut(p,v,...) for the v.x()<0 case
36 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
37 // 30.11.04 V.Grichine: modifications in SurfaceNormal for edges/vertices and
38 // in constructor with vertices
39 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
40 // 18.11.99 V.Grichine: kUndef was added to ESide
41 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
42 // 21.03.95 P.Kent: Modified for `tolerant' geom
43 //
45 
46 #include "G4Para.hh"
47 
48 #include "G4VoxelLimits.hh"
49 #include "G4AffineTransform.hh"
50 #include "Randomize.hh"
51 
52 #include "G4VPVParameterisation.hh"
53 
54 #include "G4VGraphicsScene.hh"
55 
56 using namespace CLHEP;
57 
58 // Private enum: Not for external use
59 
60 enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
61 
62 // used internally for normal routine
63 
64 enum ENSide {kNZ,kNX,kNY};
65 
67 //
68 // Constructor - check and set half-widths
69 
70 void G4Para::SetAllParameters( G4double pDx, G4double pDy, G4double pDz,
71  G4double pAlpha, G4double pTheta, G4double pPhi )
72 {
73  if ( pDx > 0 && pDy > 0 && pDz > 0 )
74  {
75  fDx = pDx;
76  fDy = pDy;
77  fDz = pDz;
78  fTalpha = std::tan(pAlpha);
79  fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
80  fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
81  }
82  else
83  {
84  std::ostringstream message;
85  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
86  << " pDx, pDy, pDz = "
87  << pDx << ", " << pDy << ", " << pDz;
88  G4Exception("G4Para::SetAllParameters()", "GeomSolids0002",
89  FatalException, message);
90  }
91  fCubicVolume = 0.;
92  fSurfaceArea = 0.;
93 }
94 
96 //
97 
98 G4Para::G4Para(const G4String& pName,
99  G4double pDx, G4double pDy, G4double pDz,
100  G4double pAlpha, G4double pTheta, G4double pPhi)
101  : G4CSGSolid(pName)
102 {
103  if ((pDx<=0) || (pDy<=0) || (pDz<=0))
104  {
105  std::ostringstream message;
106  message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
107  << " pDx, pDy, pDz = "
108  << pDx << ", " << pDy << ", " << pDz;
109  G4Exception("G4Para::G4Para()", "GeomSolids0002",
110  FatalException, message);
111  }
112  SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
113 }
114 
116 //
117 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
118 // which are its vertices. Checking of planarity with preparation of
119 // fPlanes[] and than calculation of other members
120 
121 G4Para::G4Para( const G4String& pName,
122  const G4ThreeVector pt[8] )
123  : G4CSGSolid(pName)
124 {
125  if (!( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() &&
126  pt[0].z()==pt[3].z() && pt[4].z()>0 && pt[4].z()==pt[5].z() &&
127  pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z() &&
128  (pt[0].z()+pt[4].z())==0 &&
129  pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y() &&
130  pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y() &&
131  ( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) == 0 &&
132  ( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() ) == 0) )
133  {
134  std::ostringstream message;
135  message << "Invalid vertice coordinates for Solid: " << GetName();
136  G4Exception("G4Para::G4Para()", "GeomSolids0002",
137  FatalException, message);
138  }
139  fDx = ((pt[3]).x()-(pt[2]).x())*0.5;
140  fDy = ((pt[2]).y()-(pt[1]).y())*0.5;
141  fDz = (pt[7]).z();
142  fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
143  fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
144  fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
145  fCubicVolume = 0.;
146  fSurfaceArea = 0.;
147 }
148 
150 //
151 // Fake default constructor - sets only member data and allocates memory
152 // for usage restricted to object persistency.
153 //
154 G4Para::G4Para( __void__& a )
155  : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.),
156  fTalpha(0.), fTthetaCphi(0.), fTthetaSphi(0.)
157 {
158 }
159 
161 //
162 
163 G4Para::~G4Para()
164 {
165 }
166 
168 //
169 // Copy constructor
170 
171 G4Para::G4Para(const G4Para& rhs)
172  : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
173  fTalpha(rhs.fTalpha), fTthetaCphi(rhs.fTthetaCphi),
174  fTthetaSphi(rhs.fTthetaSphi)
175 {
176 }
177 
179 //
180 // Assignment operator
181 
182 G4Para& G4Para::operator = (const G4Para& rhs)
183 {
184  // Check assignment to self
185  //
186  if (this == &rhs) { return *this; }
187 
188  // Copy base class data
189  //
190  G4CSGSolid::operator=(rhs);
191 
192  // Copy data
193  //
194  fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz;
195  fTalpha = rhs.fTalpha; fTthetaCphi = rhs.fTthetaCphi;
196  fTthetaSphi = rhs.fTthetaSphi;
197 
198  return *this;
199 }
200 
202 //
203 // Dispatch to parameterisation for replication mechanism dimension
204 // computation & modification.
205 
206 void G4Para::ComputeDimensions( G4VPVParameterisation* p,
207  const G4int n,
208  const G4VPhysicalVolume* pRep )
209 {
210  p->ComputeDimensions(*this,n,pRep);
211 }
212 
213 
215 //
216 // Calculate extent under transform and specified limit
217 
218 G4bool G4Para::CalculateExtent( const EAxis pAxis,
219  const G4VoxelLimits& pVoxelLimit,
220  const G4AffineTransform& pTransform,
221  G4double& pMin, G4double& pMax ) const
222 {
223  G4bool flag;
224 
225  if (!pTransform.IsRotated())
226  {
227  // Special case handling for unrotated trapezoids
228  // Compute z/x/y/ mins and maxs respecting limits, with early returns
229  // if outside limits. Then switch() on pAxis
230 
231  G4int i ;
232  G4double xoffset,xMin,xMax;
233  G4double yoffset,yMin,yMax;
234  G4double zoffset,zMin,zMax;
235  G4double temp[8] ; // some points for intersection with zMin/zMax
236 
237  xoffset=pTransform.NetTranslation().x();
238  yoffset=pTransform.NetTranslation().y();
239  zoffset=pTransform.NetTranslation().z();
240 
241  G4ThreeVector pt[8]; // vertices after translation
242  pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
243  yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
244  pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
245  yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
246  pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
247  yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
248  pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
249  yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
250  pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
251  yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
252  pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
253  yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
254  pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
255  yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
256  pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
257  yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
258  zMin=zoffset-fDz;
259  zMax=zoffset+fDz;
260  if ( pVoxelLimit.IsZLimited() )
261  {
262  if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
263  || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
264  {
265  return false;
266  }
267  else
268  {
269  if (zMin<pVoxelLimit.GetMinZExtent())
270  {
271  zMin=pVoxelLimit.GetMinZExtent();
272  }
273  if (zMax>pVoxelLimit.GetMaxZExtent())
274  {
275  zMax=pVoxelLimit.GetMaxZExtent();
276  }
277  }
278  }
279 
280  temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())
281  *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
282  temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())
283  *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
284  temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())
285  *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
286  temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())
287  *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
288  yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy - fDy ;
289  yMin = -yMax ;
290  for(i=0;i<4;i++)
291  {
292  if(temp[i] > yMax) yMax = temp[i] ;
293  if(temp[i] < yMin) yMin = temp[i] ;
294  }
295 
296  if (pVoxelLimit.IsYLimited())
297  {
298  if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
299  || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
300  {
301  return false;
302  }
303  else
304  {
305  if (yMin<pVoxelLimit.GetMinYExtent())
306  {
307  yMin=pVoxelLimit.GetMinYExtent();
308  }
309  if (yMax>pVoxelLimit.GetMaxYExtent())
310  {
311  yMax=pVoxelLimit.GetMaxYExtent();
312  }
313  }
314  }
315 
316  temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
317  *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
318  temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
319  *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
320  temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
321  *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
322  temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
323  *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
324  temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
325  *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
326  temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
327  *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
328  temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
329  *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
330  temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
331  *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
332 
333  xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
334  xMin = -xMax ;
335  for(i=0;i<8;i++)
336  {
337  if(temp[i] > xMax) xMax = temp[i] ;
338  if(temp[i] < xMin) xMin = temp[i] ;
339  }
340  // xMax/Min = f(yMax/Min) ?
341  if (pVoxelLimit.IsXLimited())
342  {
343  if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
344  || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
345  {
346  return false;
347  }
348  else
349  {
350  if (xMin<pVoxelLimit.GetMinXExtent())
351  {
352  xMin=pVoxelLimit.GetMinXExtent();
353  }
354  if (xMax>pVoxelLimit.GetMaxXExtent())
355  {
356  xMax=pVoxelLimit.GetMaxXExtent();
357  }
358  }
359  }
360 
361  switch (pAxis)
362  {
363  case kXAxis:
364  pMin=xMin;
365  pMax=xMax;
366  break;
367  case kYAxis:
368  pMin=yMin;
369  pMax=yMax;
370  break;
371  case kZAxis:
372  pMin=zMin;
373  pMax=zMax;
374  break;
375  default:
376  break;
377  }
378 
379  pMin-=kCarTolerance;
380  pMax+=kCarTolerance;
381  flag = true;
382  }
383  else
384  {
385  // General rotated case - create and clip mesh to boundaries
386 
387  G4bool existsAfterClip=false;
388  G4ThreeVectorList *vertices;
389 
390  pMin=+kInfinity;
391  pMax=-kInfinity;
392 
393  // Calculate rotated vertex coordinates
394 
395  vertices=CreateRotatedVertices(pTransform);
396  ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
397  ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
398  ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
399 
400  if (pMin!=kInfinity||pMax!=-kInfinity)
401  {
402  existsAfterClip=true;
403 
404  // Add 2*tolerance to avoid precision troubles
405  //
406  pMin-=kCarTolerance;
407  pMax+=kCarTolerance;
408  }
409  else
410  {
411  // Check for case where completely enveloping clipping volume
412  // If point inside then we are confident that the solid completely
413  // envelopes the clipping volume. Hence set min/max extents according
414  // to clipping volume extents along the specified axis.
415 
416  G4ThreeVector clipCentre(
417  (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
418  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
419  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
420 
421  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
422  {
423  existsAfterClip=true;
424  pMin=pVoxelLimit.GetMinExtent(pAxis);
425  pMax=pVoxelLimit.GetMaxExtent(pAxis);
426  }
427  }
428  delete vertices ; // 'new' in the function called
429  flag = existsAfterClip ;
430  }
431  return flag;
432 }
433 
435 //
436 // Check in p is inside/on surface/outside solid
437 
438 EInside G4Para::Inside( const G4ThreeVector& p ) const
439 {
440  G4double xt, yt, yt1;
441  EInside in = kOutside;
442 
443  yt1 = p.y() - fTthetaSphi*p.z();
444  yt = std::fabs(yt1) ;
445 
446  // xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
447 
448  xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
449 
450  if ( std::fabs( p.z() ) <= fDz - kCarTolerance*0.5)
451  {
452  if (yt <= fDy - kCarTolerance*0.5)
453  {
454  if ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
455  else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
456  }
457  else if ( yt <= fDy + kCarTolerance*0.5)
458  {
459  if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
460  }
461  }
462  else if ( std::fabs(p.z()) <= fDz + kCarTolerance*0.5 )
463  {
464  if ( yt <= fDy + kCarTolerance*0.5)
465  {
466  if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
467  }
468  }
469  return in;
470 }
471 
473 //
474 // Calculate side nearest to p, and return normal
475 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
476 
477 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
478 {
479  G4ThreeVector norm, sumnorm(0.,0.,0.);
480  G4int noSurfaces = 0;
481  G4double distx,disty,distz;
482  G4double newpx,newpy,xshift;
483  G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
484  G4double tntheta,cosntheta; // tan and cos of normal's theta component
485  G4double ycomp;
486  G4double delta = 0.5*kCarTolerance;
487 
488  newpx = p.x()-fTthetaCphi*p.z();
489  newpy = p.y()-fTthetaSphi*p.z();
490 
491  calpha = 1/std::sqrt(1+fTalpha*fTalpha);
492  salpha = -calpha*fTalpha; // NOTE: using MINUS std::sin(alpha)
493 
494  // xshift = newpx*calpha+newpy*salpha;
495  xshift = newpx - newpy*fTalpha;
496 
497  // distx = std::fabs(std::fabs(xshift)-fDx*calpha);
498  distx = std::fabs(std::fabs(xshift)-fDx);
499  disty = std::fabs(std::fabs(newpy)-fDy);
500  distz = std::fabs(std::fabs(p.z())-fDz);
501 
502  tntheta = fTthetaCphi*calpha + fTthetaSphi*salpha;
503  cosntheta = 1/std::sqrt(1+tntheta*tntheta);
504  ycomp = 1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
505 
506  G4ThreeVector nX = G4ThreeVector( calpha*cosntheta,
507  salpha*cosntheta,
508  -tntheta*cosntheta);
509  G4ThreeVector nY = G4ThreeVector( 0, ycomp,-fTthetaSphi*ycomp);
510  G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
511 
512  if (distx <= delta)
513  {
514  noSurfaces ++;
515  if ( xshift >= 0.) {sumnorm += nX;}
516  else {sumnorm -= nX;}
517  }
518  if (disty <= delta)
519  {
520  noSurfaces ++;
521  if ( newpy >= 0.) {sumnorm += nY;}
522  else {sumnorm -= nY;}
523  }
524  if (distz <= delta)
525  {
526  noSurfaces ++;
527  if ( p.z() >= 0.) {sumnorm += nZ;}
528  else {sumnorm -= nZ;}
529  }
530  if ( noSurfaces == 0 )
531  {
532 #ifdef G4CSGDEBUG
533  G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
534  JustWarning, "Point p is not on surface !?" );
535 #endif
536  norm = ApproxSurfaceNormal(p);
537  }
538  else if ( noSurfaces == 1 ) {norm = sumnorm;}
539  else {norm = sumnorm.unit();}
540 
541  return norm;
542 }
543 
544 
546 //
547 // Algorithm for SurfaceNormal() following the original specification
548 // for points not on the surface
549 
550 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
551 {
552  ENSide side;
553  G4ThreeVector norm;
554  G4double distx,disty,distz;
555  G4double newpx,newpy,xshift;
556  G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
557  G4double tntheta,cosntheta; // tan and cos of normal's theta component
558  G4double ycomp;
559 
560  newpx=p.x()-fTthetaCphi*p.z();
561  newpy=p.y()-fTthetaSphi*p.z();
562 
563  calpha=1/std::sqrt(1+fTalpha*fTalpha);
564  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
565 
566  xshift=newpx*calpha+newpy*salpha;
567 
568  distx=std::fabs(std::fabs(xshift)-fDx*calpha);
569  disty=std::fabs(std::fabs(newpy)-fDy);
570  distz=std::fabs(std::fabs(p.z())-fDz);
571 
572  if (distx<disty)
573  {
574  if (distx<distz) {side=kNX;}
575  else {side=kNZ;}
576  }
577  else
578  {
579  if (disty<distz) {side=kNY;}
580  else {side=kNZ;}
581  }
582 
583  switch (side)
584  {
585  case kNX:
586  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
587  if (xshift<0)
588  {
589  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
590  }
591  else
592  {
593  cosntheta=1/std::sqrt(1+tntheta*tntheta);
594  }
595  norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
596  break;
597  case kNY:
598  if (newpy<0)
599  {
600  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
601  }
602  else
603  {
604  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
605  }
606  norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
607  break;
608  case kNZ: // Closest to Z
609  if (p.z()>=0)
610  {
611  norm=G4ThreeVector(0,0,1);
612  }
613  else
614  {
615  norm=G4ThreeVector(0,0,-1);
616  }
617  break;
618  }
619  return norm;
620 }
621 
623 //
624 // Calculate distance to shape from outside
625 // - return kInfinity if no intersection
626 //
627 // ALGORITHM:
628 // For each component, calculate pair of minimum and maximum intersection
629 // values for which the particle is in the extent of the shape
630 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
631 // - Z plane intersectin uses tolerance
632 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
633 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
634 // cases)
635 // - Note: XZ and YZ planes each divide space into four regions,
636 // characterised by ss1 ss2
637 
638 G4double G4Para::DistanceToIn( const G4ThreeVector& p,
639  const G4ThreeVector& v ) const
640 {
641  G4double snxt; // snxt = default return value
642  G4double smin,smax;
643  G4double tmin,tmax;
644  G4double yt,vy,xt,vx;
645  G4double max;
646  //
647  // Z Intersection range
648  //
649  if (v.z()>0)
650  {
651  max=fDz-p.z();
652  if (max>kCarTolerance*0.5)
653  {
654  smax=max/v.z();
655  smin=(-fDz-p.z())/v.z();
656  }
657  else
658  {
659  return snxt=kInfinity;
660  }
661  }
662  else if (v.z()<0)
663  {
664  max=-fDz-p.z();
665  if (max<-kCarTolerance*0.5)
666  {
667  smax=max/v.z();
668  smin=(fDz-p.z())/v.z();
669  }
670  else
671  {
672  return snxt=kInfinity;
673  }
674  }
675  else
676  {
677  if (std::fabs(p.z())<=fDz) // Inside
678  {
679  smin=0;
680  smax=kInfinity;
681  }
682  else
683  {
684  return snxt=kInfinity;
685  }
686  }
687 
688  //
689  // Y G4Parallel planes intersection
690  //
691 
692  yt=p.y()-fTthetaSphi*p.z();
693  vy=v.y()-fTthetaSphi*v.z();
694 
695  if (vy>0)
696  {
697  max=fDy-yt;
698  if (max>kCarTolerance*0.5)
699  {
700  tmax=max/vy;
701  tmin=(-fDy-yt)/vy;
702  }
703  else
704  {
705  return snxt=kInfinity;
706  }
707  }
708  else if (vy<0)
709  {
710  max=-fDy-yt;
711  if (max<-kCarTolerance*0.5)
712  {
713  tmax=max/vy;
714  tmin=(fDy-yt)/vy;
715  }
716  else
717  {
718  return snxt=kInfinity;
719  }
720  }
721  else
722  {
723  if (std::fabs(yt)<=fDy)
724  {
725  tmin=0;
726  tmax=kInfinity;
727  }
728  else
729  {
730  return snxt=kInfinity;
731  }
732  }
733 
734  // Re-Calc valid intersection range
735  //
736  if (tmin>smin) smin=tmin;
737  if (tmax<smax) smax=tmax;
738  if (smax<=smin)
739  {
740  return snxt=kInfinity;
741  }
742  else
743  {
744  //
745  // X G4Parallel planes intersection
746  //
747  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
748  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
749  if (vx>0)
750  {
751  max=fDx-xt;
752  if (max>kCarTolerance*0.5)
753  {
754  tmax=max/vx;
755  tmin=(-fDx-xt)/vx;
756  }
757  else
758  {
759  return snxt=kInfinity;
760  }
761  }
762  else if (vx<0)
763  {
764  max=-fDx-xt;
765  if (max<-kCarTolerance*0.5)
766  {
767  tmax=max/vx;
768  tmin=(fDx-xt)/vx;
769  }
770  else
771  {
772  return snxt=kInfinity;
773  }
774  }
775  else
776  {
777  if (std::fabs(xt)<=fDx)
778  {
779  tmin=0;
780  tmax=kInfinity;
781  }
782  else
783  {
784  return snxt=kInfinity;
785  }
786  }
787  if (tmin>smin) smin=tmin;
788  if (tmax<smax) smax=tmax;
789  }
790 
791  if (smax>0&&smin<smax)
792  {
793  if (smin>0)
794  {
795  snxt=smin;
796  }
797  else
798  {
799  snxt=0;
800  }
801  }
802  else
803  {
804  snxt=kInfinity;
805  }
806  return snxt;
807 }
808 
810 //
811 // Calculate exact shortest distance to any boundary from outside
812 // - Returns 0 is point inside
813 
814 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
815 {
816  G4double safe=0.0;
817  G4double distz1,distz2,disty1,disty2,distx1,distx2;
818  G4double trany,cosy,tranx,cosx;
819 
820  // Z planes
821  //
822  distz1=p.z()-fDz;
823  distz2=-fDz-p.z();
824  if (distz1>distz2)
825  {
826  safe=distz1;
827  }
828  else
829  {
830  safe=distz2;
831  }
832 
833  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
834 
835  // Transformed x into `box' system
836  //
837  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
838  disty1=(trany-fDy)*cosy;
839  disty2=(-fDy-trany)*cosy;
840 
841  if (disty1>safe) safe=disty1;
842  if (disty2>safe) safe=disty2;
843 
844  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
845  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
846  distx1=(tranx-fDx)*cosx;
847  distx2=(-fDx-tranx)*cosx;
848 
849  if (distx1>safe) safe=distx1;
850  if (distx2>safe) safe=distx2;
851 
852  if (safe<0) safe=0;
853  return safe;
854 }
855 
857 //
858 // Calculate distance to surface of shape from inside
859 // Calculate distance to x/y/z planes - smallest is exiting distance
860 
861 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
862  const G4bool calcNorm,
863  G4bool *validNorm, G4ThreeVector *n) const
864 {
865  ESide side = kUndef;
866  G4double snxt; // snxt = return value
867  G4double max,tmax;
868  G4double yt,vy,xt,vx;
869 
870  G4double ycomp,calpha,salpha,tntheta,cosntheta;
871 
872  //
873  // Z Intersections
874  //
875 
876  if (v.z()>0)
877  {
878  max=fDz-p.z();
879  if (max>kCarTolerance*0.5)
880  {
881  snxt=max/v.z();
882  side=kPZ;
883  }
884  else
885  {
886  if (calcNorm)
887  {
888  *validNorm=true;
889  *n=G4ThreeVector(0,0,1);
890  }
891  return snxt=0;
892  }
893  }
894  else if (v.z()<0)
895  {
896  max=-fDz-p.z();
897  if (max<-kCarTolerance*0.5)
898  {
899  snxt=max/v.z();
900  side=kMZ;
901  }
902  else
903  {
904  if (calcNorm)
905  {
906  *validNorm=true;
907  *n=G4ThreeVector(0,0,-1);
908  }
909  return snxt=0;
910  }
911  }
912  else
913  {
914  snxt=kInfinity;
915  }
916 
917  //
918  // Y plane intersection
919  //
920 
921  yt=p.y()-fTthetaSphi*p.z();
922  vy=v.y()-fTthetaSphi*v.z();
923 
924  if (vy>0)
925  {
926  max=fDy-yt;
927  if (max>kCarTolerance*0.5)
928  {
929  tmax=max/vy;
930  if (tmax<snxt)
931  {
932  snxt=tmax;
933  side=kPY;
934  }
935  }
936  else
937  {
938  if (calcNorm)
939  {
940  *validNorm=true; // Leaving via plus Y
941  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
942  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
943  }
944  return snxt=0;
945  }
946  }
947  else if (vy<0)
948  {
949  max=-fDy-yt;
950  if (max<-kCarTolerance*0.5)
951  {
952  tmax=max/vy;
953  if (tmax<snxt)
954  {
955  snxt=tmax;
956  side=kMY;
957  }
958  }
959  else
960  {
961  if (calcNorm)
962  {
963  *validNorm=true; // Leaving via minus Y
964  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
965  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
966  }
967  return snxt=0;
968  }
969  }
970 
971  //
972  // X plane intersection
973  //
974 
975  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
976  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
977  if (vx>0)
978  {
979  max=fDx-xt;
980  if (max>kCarTolerance*0.5)
981  {
982  tmax=max/vx;
983  if (tmax<snxt)
984  {
985  snxt=tmax;
986  side=kPX;
987  }
988  }
989  else
990  {
991  if (calcNorm)
992  {
993  *validNorm=true; // Leaving via plus X
994  calpha=1/std::sqrt(1+fTalpha*fTalpha);
995  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
996  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
997  cosntheta=1/std::sqrt(1+tntheta*tntheta);
998  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
999  }
1000  return snxt=0;
1001  }
1002  }
1003  else if (vx<0)
1004  {
1005  max=-fDx-xt;
1006  if (max<-kCarTolerance*0.5)
1007  {
1008  tmax=max/vx;
1009  if (tmax<snxt)
1010  {
1011  snxt=tmax;
1012  side=kMX;
1013  }
1014  }
1015  else
1016  {
1017  if (calcNorm)
1018  {
1019  *validNorm=true; // Leaving via minus X
1020  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1021  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1022  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1023  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
1024  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1025  }
1026  return snxt=0;
1027  }
1028  }
1029 
1030  if (calcNorm)
1031  {
1032  *validNorm=true;
1033  switch (side)
1034  {
1035  case kMZ:
1036  *n=G4ThreeVector(0,0,-1);
1037  break;
1038  case kPZ:
1039  *n=G4ThreeVector(0,0,1);
1040  break;
1041  case kMY:
1042  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
1043  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
1044  break;
1045  case kPY:
1046  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
1047  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
1048  break;
1049  case kMX:
1050  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1051  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1052  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1053  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
1054  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1055  break;
1056  case kPX:
1057  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1058  salpha=-calpha*fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1059  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1060  cosntheta=1/std::sqrt(1+tntheta*tntheta);
1061  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1062  break;
1063  default:
1064  DumpInfo();
1065  G4Exception("G4Para::DistanceToOut(p,v,..)",
1066  "GeomSolids1002",JustWarning,
1067  "Undefined side for valid surface normal to solid.");
1068  break;
1069  }
1070  }
1071  return snxt;
1072 }
1073 
1075 //
1076 // Calculate exact shortest distance to any boundary from inside
1077 // - Returns 0 is point outside
1078 
1079 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
1080 {
1081  G4double safe=0.0;
1082  G4double distz1,distz2,disty1,disty2,distx1,distx2;
1083  G4double trany,cosy,tranx,cosx;
1084 
1085 #ifdef G4CSGDEBUG
1086  if( Inside(p) == kOutside )
1087  {
1088  G4int oldprc = G4cout.precision(16) ;
1089  G4cout << G4endl ;
1090  DumpInfo();
1091  G4cout << "Position:" << G4endl << G4endl ;
1092  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1093  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1094  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1095  G4cout.precision(oldprc) ;
1096  G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
1097  JustWarning, "Point p is outside !?" );
1098  }
1099 #endif
1100 
1101  // Z planes
1102  //
1103  distz1=fDz-p.z();
1104  distz2=fDz+p.z();
1105  if (distz1<distz2)
1106  {
1107  safe=distz1;
1108  }
1109  else
1110  {
1111  safe=distz2;
1112  }
1113 
1114  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
1115 
1116  // Transformed x into `box' system
1117  //
1118  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
1119  disty1=(fDy-trany)*cosy;
1120  disty2=(fDy+trany)*cosy;
1121 
1122  if (disty1<safe) safe=disty1;
1123  if (disty2<safe) safe=disty2;
1124 
1125  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
1126  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
1127  distx1=(fDx-tranx)*cosx;
1128  distx2=(fDx+tranx)*cosx;
1129 
1130  if (distx1<safe) safe=distx1;
1131  if (distx2<safe) safe=distx2;
1132 
1133  if (safe<0) safe=0;
1134  return safe;
1135 }
1136 
1138 //
1139 // Create a List containing the transformed vertices
1140 // Ordering [0-3] -fDz cross section
1141 // [4-7] +fDz cross section such that [0] is below [4],
1142 // [1] below [5] etc.
1143 // Note:
1144 // Caller has deletion resposibility
1145 
1146 G4ThreeVectorList*
1147 G4Para::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1148 {
1149  G4ThreeVectorList *vertices;
1150  vertices=new G4ThreeVectorList();
1151  if (vertices)
1152  {
1153  vertices->reserve(8);
1154  G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
1155  -fDz*fTthetaSphi-fDy, -fDz);
1156  G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
1157  -fDz*fTthetaSphi-fDy, -fDz);
1158  G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
1159  -fDz*fTthetaSphi+fDy, -fDz);
1160  G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
1161  -fDz*fTthetaSphi+fDy, -fDz);
1162  G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
1163  +fDz*fTthetaSphi-fDy, +fDz);
1164  G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
1165  +fDz*fTthetaSphi-fDy, +fDz);
1166  G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
1167  +fDz*fTthetaSphi+fDy, +fDz);
1168  G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
1169  +fDz*fTthetaSphi+fDy, +fDz);
1170 
1171  vertices->push_back(pTransform.TransformPoint(vertex0));
1172  vertices->push_back(pTransform.TransformPoint(vertex1));
1173  vertices->push_back(pTransform.TransformPoint(vertex2));
1174  vertices->push_back(pTransform.TransformPoint(vertex3));
1175  vertices->push_back(pTransform.TransformPoint(vertex4));
1176  vertices->push_back(pTransform.TransformPoint(vertex5));
1177  vertices->push_back(pTransform.TransformPoint(vertex6));
1178  vertices->push_back(pTransform.TransformPoint(vertex7));
1179  }
1180  else
1181  {
1182  DumpInfo();
1183  G4Exception("G4Para::CreateRotatedVertices()",
1184  "GeomSolids0003", FatalException,
1185  "Error in allocation of vertices. Out of memory !");
1186  }
1187  return vertices;
1188 }
1189 
1191 //
1192 // GetEntityType
1193 
1194 G4GeometryType G4Para::GetEntityType() const
1195 {
1196  return G4String("G4Para");
1197 }
1198 
1200 //
1201 // Make a clone of the object
1202 //
1203 G4VSolid* G4Para::Clone() const
1204 {
1205  return new G4Para(*this);
1206 }
1207 
1209 //
1210 // Stream object contents to an output stream
1211 
1212 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
1213 {
1214  G4int oldprc = os.precision(16);
1215  os << "-----------------------------------------------------------\n"
1216  << " *** Dump for solid - " << GetName() << " ***\n"
1217  << " ===================================================\n"
1218  << " Solid type: G4Para\n"
1219  << " Parameters: \n"
1220  << " half length X: " << fDx/mm << " mm \n"
1221  << " half length Y: " << fDy/mm << " mm \n"
1222  << " half length Z: " << fDz/mm << " mm \n"
1223  << " std::tan(alpha) : " << fTalpha/degree << " degrees \n"
1224  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree
1225  << " degrees \n"
1226  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree
1227  << " degrees \n"
1228  << "-----------------------------------------------------------\n";
1229  os.precision(oldprc);
1230 
1231  return os;
1232 }
1233 
1235 //
1236 // GetPointOnPlane
1237 // Auxiliary method for Get Point on Surface
1238 //
1239 
1240 G4ThreeVector G4Para::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1241  G4ThreeVector p2, G4ThreeVector p3,
1242  G4double& area) const
1243 {
1244  G4double lambda1, lambda2, chose, aOne, aTwo;
1245  G4ThreeVector t, u, v, w, Area, normal;
1246 
1247  t = p1 - p0;
1248  u = p2 - p1;
1249  v = p3 - p2;
1250  w = p0 - p3;
1251 
1252  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1253  w.z()*v.x() - w.x()*v.z(),
1254  w.x()*v.y() - w.y()*v.x());
1255 
1256  aOne = 0.5*Area.mag();
1257 
1258  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1259  t.z()*u.x() - t.x()*u.z(),
1260  t.x()*u.y() - t.y()*u.x());
1261 
1262  aTwo = 0.5*Area.mag();
1263 
1264  area = aOne + aTwo;
1265 
1266  chose = RandFlat::shoot(0.,aOne+aTwo);
1267 
1268  if( (chose>=0.) && (chose < aOne) )
1269  {
1270  lambda1 = RandFlat::shoot(0.,1.);
1271  lambda2 = RandFlat::shoot(0.,lambda1);
1272  return (p2+lambda1*v+lambda2*w);
1273  }
1274 
1275  // else
1276 
1277  lambda1 = RandFlat::shoot(0.,1.);
1278  lambda2 = RandFlat::shoot(0.,lambda1);
1279  return (p0+lambda1*t+lambda2*u);
1280 }
1281 
1283 //
1284 // GetPointOnSurface
1285 //
1286 // Return a point (G4ThreeVector) randomly and uniformly
1287 // selected on the solid surface
1288 
1289 G4ThreeVector G4Para::GetPointOnSurface() const
1290 {
1291  G4ThreeVector One, Two, Three, Four, Five, Six;
1292  G4ThreeVector pt[8] ;
1293  G4double chose, aOne, aTwo, aThree, aFour, aFive, aSix;
1294 
1295  pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
1296  -fDz*fTthetaSphi-fDy, -fDz);
1297  pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
1298  -fDz*fTthetaSphi-fDy, -fDz);
1299  pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
1300  -fDz*fTthetaSphi+fDy, -fDz);
1301  pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
1302  -fDz*fTthetaSphi+fDy, -fDz);
1303  pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
1304  +fDz*fTthetaSphi-fDy, +fDz);
1305  pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
1306  +fDz*fTthetaSphi-fDy, +fDz);
1307  pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
1308  +fDz*fTthetaSphi+fDy, +fDz);
1309  pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
1310  +fDz*fTthetaSphi+fDy, +fDz);
1311 
1312  // make sure we provide the points in a clockwise fashion
1313 
1314  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1315  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1316  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1317  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1318  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1319  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1320 
1321  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1322 
1323  if( (chose>=0.) && (chose<aOne) )
1324  { return One; }
1325  else if(chose>=aOne && chose<aOne+aTwo)
1326  { return Two; }
1327  else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
1328  { return Three; }
1329  else if(chose>=aOne+aTwo+aThree && chose<aOne+aTwo+aThree+aFour)
1330  { return Four; }
1331  else if(chose>=aOne+aTwo+aThree+aFour && chose<aOne+aTwo+aThree+aFour+aFive)
1332  { return Five; }
1333  return Six;
1334 }
1335 
1337 //
1338 // Methods for visualisation
1339 
1340 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1341 {
1342  scene.AddSolid (*this);
1343 }
1344 
1345 G4Polyhedron* G4Para::CreatePolyhedron () const
1346 {
1347  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1348  G4double alpha = std::atan(fTalpha);
1349  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1350  +fTthetaSphi*fTthetaSphi));
1351 
1352  return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
1353 }
G4VSolid::GetName
G4String GetName() const
G4Para::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Para.cc:861
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
kMX
Definition: G4Para.cc:60
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:68
G4Para::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Para.cc:1289
G4Para
Definition: G4Para.hh:77
G4Para::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Para.cc:206
G4Para::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Para.cc:438
G4Para::SetAllParameters
void SetAllParameters(G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:70
G4Para::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Para.cc:218
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VGraphicsScene.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4VoxelLimits.hh
G4Para::fTalpha
G4double fTalpha
Definition: G4Para.hh:180
G4Para::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Para.cc:550
z
G4double z
Definition: TRTMaterials.hh:39
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4Para::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Para.cc:477
G4Para::fDx
G4double fDx
Definition: G4Para.hh:179
kMY
Definition: G4Para.cc:60
G4Para::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Para.cc:638
G4AffineTransform::IsRotated
G4bool IsRotated() const
w
const G4double w[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:44
G4Para::G4Para
G4Para(const G4String &pName, G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:98
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
kNZ
Definition: G4Para.cc:64
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4Para::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Para.cc:1212
a
G4double a
Definition: TRTMaterials.hh:39
G4Para::operator=
G4Para & operator=(const G4Para &rhs)
Definition: G4Para.cc:182
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
smax
const G4int smax
Definition: G4UrQMD1_3Interface.hh:71
G4Para::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Para.cc:1345
G4VSolid::DumpInfo
void DumpInfo() const
G4Para::fDz
G4double fDz
Definition: G4Para.hh:179
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
kPZ
Definition: G4Para.cc:60
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
kNY
Definition: G4Para.cc:64
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4cout
G4GLOB_DLL std::ostream G4cout
kNX
Definition: G4Para.cc:64
kPX
Definition: G4Para.cc:60
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VPVParameterisation.hh
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
kUndef
Definition: G4Para.cc:60
Randomize.hh
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4Para::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Para.cc:1194
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
ENSide
ENSide
Definition: G4Para.cc:64
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4Para::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Para.cc:1340
kMZ
Definition: G4Para.cc:60
G4AffineTransform.hh
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
EInside
EInside
Definition: geomdefs.hh:58
G4Para::fTthetaCphi
G4double fTthetaCphi
Definition: G4Para.hh:180
EAxis
EAxis
Definition: geomdefs.hh:54
G4Para::fTthetaSphi
G4double fTthetaSphi
Definition: G4Para.hh:180
x
const G4double x[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:56
FatalException
Definition: G4ExceptionSeverity.hh:60
G4Para::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Para.cc:1147
degree
static const double degree
Definition: G4SIunits.hh:143
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4Para::Clone
G4VSolid * Clone() const
Definition: G4Para.cc:1203
G4Para.hh
ESide
ESide
Definition: G4Cons.cc:71
G4Para::~G4Para
virtual ~G4Para()
Definition: G4Para.cc:163
G4double
double G4double
Definition: G4Types.hh:76
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:79
kPY
Definition: G4Para.cc:60
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
mm
static const double mm
Definition: G4SIunits.hh:114
G4VSolid::ClipBetweenSections
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
G4Para::fDy
G4double fDy
Definition: G4Para.hh:179
G4Para::GetPointOnPlane
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double &area) const
Definition: G4Para.cc:1240