Geant4  10.00.p03
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 84622 2014-10-17 09:43:17Z 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  if (fTalpha) {salpha = -calpha*fTalpha;} // NOTE: using MINUS std::sin(alpha)
493  else {salpha = 0.;}
494 
495  // xshift = newpx*calpha+newpy*salpha;
496  xshift = newpx - newpy*fTalpha;
497 
498  // distx = std::fabs(std::fabs(xshift)-fDx*calpha);
499  distx = std::fabs(std::fabs(xshift)-fDx);
500  disty = std::fabs(std::fabs(newpy)-fDy);
501  distz = std::fabs(std::fabs(p.z())-fDz);
502 
503  tntheta = fTthetaCphi*calpha + fTthetaSphi*salpha;
504  cosntheta = 1/std::sqrt(1+tntheta*tntheta);
505  ycomp = 1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
506 
507  G4ThreeVector nX = G4ThreeVector( calpha*cosntheta,
508  salpha*cosntheta,
509  -tntheta*cosntheta);
510  G4ThreeVector nY = G4ThreeVector( 0, ycomp,-fTthetaSphi*ycomp);
511  G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
512 
513  if (distx <= delta)
514  {
515  noSurfaces ++;
516  if ( xshift >= 0.) {sumnorm += nX;}
517  else {sumnorm -= nX;}
518  }
519  if (disty <= delta)
520  {
521  noSurfaces ++;
522  if ( newpy >= 0.) {sumnorm += nY;}
523  else {sumnorm -= nY;}
524  }
525  if (distz <= delta)
526  {
527  noSurfaces ++;
528  if ( p.z() >= 0.) {sumnorm += nZ;}
529  else {sumnorm -= nZ;}
530  }
531  if ( noSurfaces == 0 )
532  {
533 #ifdef G4CSGDEBUG
534  G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
535  JustWarning, "Point p is not on surface !?" );
536 #endif
537  norm = ApproxSurfaceNormal(p);
538  }
539  else if ( noSurfaces == 1 ) {norm = sumnorm;}
540  else {norm = sumnorm.unit();}
541 
542  return norm;
543 }
544 
545 
547 //
548 // Algorithm for SurfaceNormal() following the original specification
549 // for points not on the surface
550 
551 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
552 {
553  ENSide side;
554  G4ThreeVector norm;
555  G4double distx,disty,distz;
556  G4double newpx,newpy,xshift;
557  G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
558  G4double tntheta,cosntheta; // tan and cos of normal's theta component
559  G4double ycomp;
560 
561  newpx=p.x()-fTthetaCphi*p.z();
562  newpy=p.y()-fTthetaSphi*p.z();
563 
564  calpha=1/std::sqrt(1+fTalpha*fTalpha);
565  if (fTalpha)
566  {
567  salpha=-calpha/fTalpha; // NOTE: actually use MINUS std::sin(alpha)
568  }
569  else
570  {
571  salpha=0;
572  }
573 
574  xshift=newpx*calpha+newpy*salpha;
575 
576  distx=std::fabs(std::fabs(xshift)-fDx*calpha);
577  disty=std::fabs(std::fabs(newpy)-fDy);
578  distz=std::fabs(std::fabs(p.z())-fDz);
579 
580  if (distx<disty)
581  {
582  if (distx<distz) {side=kNX;}
583  else {side=kNZ;}
584  }
585  else
586  {
587  if (disty<distz) {side=kNY;}
588  else {side=kNZ;}
589  }
590 
591  switch (side)
592  {
593  case kNX:
594  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
595  if (xshift<0)
596  {
597  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
598  }
599  else
600  {
601  cosntheta=1/std::sqrt(1+tntheta*tntheta);
602  }
603  norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
604  break;
605  case kNY:
606  if (newpy<0)
607  {
608  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
609  }
610  else
611  {
612  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
613  }
614  norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
615  break;
616  case kNZ: // Closest to Z
617  if (p.z()>=0)
618  {
619  norm=G4ThreeVector(0,0,1);
620  }
621  else
622  {
623  norm=G4ThreeVector(0,0,-1);
624  }
625  break;
626  }
627  return norm;
628 }
629 
631 //
632 // Calculate distance to shape from outside
633 // - return kInfinity if no intersection
634 //
635 // ALGORITHM:
636 // For each component, calculate pair of minimum and maximum intersection
637 // values for which the particle is in the extent of the shape
638 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
639 // - Z plane intersectin uses tolerance
640 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
641 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
642 // cases)
643 // - Note: XZ and YZ planes each divide space into four regions,
644 // characterised by ss1 ss2
645 
646 G4double G4Para::DistanceToIn( const G4ThreeVector& p,
647  const G4ThreeVector& v ) const
648 {
649  G4double snxt; // snxt = default return value
650  G4double smin,smax;
651  G4double tmin,tmax;
652  G4double yt,vy,xt,vx;
653  G4double max;
654  //
655  // Z Intersection range
656  //
657  if (v.z()>0)
658  {
659  max=fDz-p.z();
660  if (max>kCarTolerance*0.5)
661  {
662  smax=max/v.z();
663  smin=(-fDz-p.z())/v.z();
664  }
665  else
666  {
667  return snxt=kInfinity;
668  }
669  }
670  else if (v.z()<0)
671  {
672  max=-fDz-p.z();
673  if (max<-kCarTolerance*0.5)
674  {
675  smax=max/v.z();
676  smin=(fDz-p.z())/v.z();
677  }
678  else
679  {
680  return snxt=kInfinity;
681  }
682  }
683  else
684  {
685  if (std::fabs(p.z())<=fDz) // Inside
686  {
687  smin=0;
688  smax=kInfinity;
689  }
690  else
691  {
692  return snxt=kInfinity;
693  }
694  }
695 
696  //
697  // Y G4Parallel planes intersection
698  //
699 
700  yt=p.y()-fTthetaSphi*p.z();
701  vy=v.y()-fTthetaSphi*v.z();
702 
703  if (vy>0)
704  {
705  max=fDy-yt;
706  if (max>kCarTolerance*0.5)
707  {
708  tmax=max/vy;
709  tmin=(-fDy-yt)/vy;
710  }
711  else
712  {
713  return snxt=kInfinity;
714  }
715  }
716  else if (vy<0)
717  {
718  max=-fDy-yt;
719  if (max<-kCarTolerance*0.5)
720  {
721  tmax=max/vy;
722  tmin=(fDy-yt)/vy;
723  }
724  else
725  {
726  return snxt=kInfinity;
727  }
728  }
729  else
730  {
731  if (std::fabs(yt)<=fDy)
732  {
733  tmin=0;
734  tmax=kInfinity;
735  }
736  else
737  {
738  return snxt=kInfinity;
739  }
740  }
741 
742  // Re-Calc valid intersection range
743  //
744  if (tmin>smin) smin=tmin;
745  if (tmax<smax) smax=tmax;
746  if (smax<=smin)
747  {
748  return snxt=kInfinity;
749  }
750  else
751  {
752  //
753  // X G4Parallel planes intersection
754  //
755  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
756  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
757  if (vx>0)
758  {
759  max=fDx-xt;
760  if (max>kCarTolerance*0.5)
761  {
762  tmax=max/vx;
763  tmin=(-fDx-xt)/vx;
764  }
765  else
766  {
767  return snxt=kInfinity;
768  }
769  }
770  else if (vx<0)
771  {
772  max=-fDx-xt;
773  if (max<-kCarTolerance*0.5)
774  {
775  tmax=max/vx;
776  tmin=(fDx-xt)/vx;
777  }
778  else
779  {
780  return snxt=kInfinity;
781  }
782  }
783  else
784  {
785  if (std::fabs(xt)<=fDx)
786  {
787  tmin=0;
788  tmax=kInfinity;
789  }
790  else
791  {
792  return snxt=kInfinity;
793  }
794  }
795  if (tmin>smin) smin=tmin;
796  if (tmax<smax) smax=tmax;
797  }
798 
799  if (smax>0&&smin<smax)
800  {
801  if (smin>0)
802  {
803  snxt=smin;
804  }
805  else
806  {
807  snxt=0;
808  }
809  }
810  else
811  {
812  snxt=kInfinity;
813  }
814  return snxt;
815 }
816 
818 //
819 // Calculate exact shortest distance to any boundary from outside
820 // - Returns 0 is point inside
821 
822 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
823 {
824  G4double safe=0.0;
825  G4double distz1,distz2,disty1,disty2,distx1,distx2;
826  G4double trany,cosy,tranx,cosx;
827 
828  // Z planes
829  //
830  distz1=p.z()-fDz;
831  distz2=-fDz-p.z();
832  if (distz1>distz2)
833  {
834  safe=distz1;
835  }
836  else
837  {
838  safe=distz2;
839  }
840 
841  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
842 
843  // Transformed x into `box' system
844  //
845  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
846  disty1=(trany-fDy)*cosy;
847  disty2=(-fDy-trany)*cosy;
848 
849  if (disty1>safe) safe=disty1;
850  if (disty2>safe) safe=disty2;
851 
852  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
853  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
854  distx1=(tranx-fDx)*cosx;
855  distx2=(-fDx-tranx)*cosx;
856 
857  if (distx1>safe) safe=distx1;
858  if (distx2>safe) safe=distx2;
859 
860  if (safe<0) safe=0;
861  return safe;
862 }
863 
865 //
866 // Calculate distance to surface of shape from inside
867 // Calculate distance to x/y/z planes - smallest is exiting distance
868 
869 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
870  const G4bool calcNorm,
871  G4bool *validNorm, G4ThreeVector *n) const
872 {
873  ESide side = kUndef;
874  G4double snxt; // snxt = return value
875  G4double max,tmax;
876  G4double yt,vy,xt,vx;
877 
878  G4double ycomp,calpha,salpha,tntheta,cosntheta;
879 
880  //
881  // Z Intersections
882  //
883 
884  if (v.z()>0)
885  {
886  max=fDz-p.z();
887  if (max>kCarTolerance*0.5)
888  {
889  snxt=max/v.z();
890  side=kPZ;
891  }
892  else
893  {
894  if (calcNorm)
895  {
896  *validNorm=true;
897  *n=G4ThreeVector(0,0,1);
898  }
899  return snxt=0;
900  }
901  }
902  else if (v.z()<0)
903  {
904  max=-fDz-p.z();
905  if (max<-kCarTolerance*0.5)
906  {
907  snxt=max/v.z();
908  side=kMZ;
909  }
910  else
911  {
912  if (calcNorm)
913  {
914  *validNorm=true;
915  *n=G4ThreeVector(0,0,-1);
916  }
917  return snxt=0;
918  }
919  }
920  else
921  {
922  snxt=kInfinity;
923  }
924 
925  //
926  // Y plane intersection
927  //
928 
929  yt=p.y()-fTthetaSphi*p.z();
930  vy=v.y()-fTthetaSphi*v.z();
931 
932  if (vy>0)
933  {
934  max=fDy-yt;
935  if (max>kCarTolerance*0.5)
936  {
937  tmax=max/vy;
938  if (tmax<snxt)
939  {
940  snxt=tmax;
941  side=kPY;
942  }
943  }
944  else
945  {
946  if (calcNorm)
947  {
948  *validNorm=true; // Leaving via plus Y
949  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
950  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
951  }
952  return snxt=0;
953  }
954  }
955  else if (vy<0)
956  {
957  max=-fDy-yt;
958  if (max<-kCarTolerance*0.5)
959  {
960  tmax=max/vy;
961  if (tmax<snxt)
962  {
963  snxt=tmax;
964  side=kMY;
965  }
966  }
967  else
968  {
969  if (calcNorm)
970  {
971  *validNorm=true; // Leaving via minus Y
972  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
973  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
974  }
975  return snxt=0;
976  }
977  }
978 
979  //
980  // X plane intersection
981  //
982 
983  xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
984  vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
985  if (vx>0)
986  {
987  max=fDx-xt;
988  if (max>kCarTolerance*0.5)
989  {
990  tmax=max/vx;
991  if (tmax<snxt)
992  {
993  snxt=tmax;
994  side=kPX;
995  }
996  }
997  else
998  {
999  if (calcNorm)
1000  {
1001  *validNorm=true; // Leaving via plus X
1002  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1003  if (fTalpha)
1004  {
1005  salpha=-calpha/fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1006  }
1007  else
1008  {
1009  salpha=0;
1010  }
1011  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1012  cosntheta=1/std::sqrt(1+tntheta*tntheta);
1013  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1014  }
1015  return snxt=0;
1016  }
1017  }
1018  else if (vx<0)
1019  {
1020  max=-fDx-xt;
1021  if (max<-kCarTolerance*0.5)
1022  {
1023  tmax=max/vx;
1024  if (tmax<snxt)
1025  {
1026  snxt=tmax;
1027  side=kMX;
1028  }
1029  }
1030  else
1031  {
1032  if (calcNorm)
1033  {
1034  *validNorm=true; // Leaving via minus X
1035  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1036  if (fTalpha)
1037  {
1038  salpha=-calpha/fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1039  }
1040  else
1041  {
1042  salpha=0;
1043  }
1044  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1045  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
1046  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1047  }
1048  return snxt=0;
1049  }
1050  }
1051 
1052  if (calcNorm)
1053  {
1054  *validNorm=true;
1055  switch (side)
1056  {
1057  case kMZ:
1058  *n=G4ThreeVector(0,0,-1);
1059  break;
1060  case kPZ:
1061  *n=G4ThreeVector(0,0,1);
1062  break;
1063  case kMY:
1064  ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
1065  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
1066  break;
1067  case kPY:
1068  ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
1069  *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
1070  break;
1071  case kMX:
1072  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1073  if (fTalpha)
1074  {
1075  salpha=-calpha/fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1076  }
1077  else
1078  {
1079  salpha=0;
1080  }
1081  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1082  cosntheta=-1/std::sqrt(1+tntheta*tntheta);
1083  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1084  break;
1085  case kPX:
1086  calpha=1/std::sqrt(1+fTalpha*fTalpha);
1087  if (fTalpha)
1088  {
1089  salpha=-calpha/fTalpha; // NOTE: actually use MINUS std::sin(alpha)
1090  }
1091  else
1092  {
1093  salpha=0;
1094  }
1095  tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
1096  cosntheta=1/std::sqrt(1+tntheta*tntheta);
1097  *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
1098  break;
1099  default:
1100  DumpInfo();
1101  G4Exception("G4Para::DistanceToOut(p,v,..)",
1102  "GeomSolids1002",JustWarning,
1103  "Undefined side for valid surface normal to solid.");
1104  break;
1105  }
1106  }
1107  return snxt;
1108 }
1109 
1111 //
1112 // Calculate exact shortest distance to any boundary from inside
1113 // - Returns 0 is point outside
1114 
1115 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
1116 {
1117  G4double safe=0.0;
1118  G4double distz1,distz2,disty1,disty2,distx1,distx2;
1119  G4double trany,cosy,tranx,cosx;
1120 
1121 #ifdef G4CSGDEBUG
1122  if( Inside(p) == kOutside )
1123  {
1124  G4int oldprc = G4cout.precision(16) ;
1125  G4cout << G4endl ;
1126  DumpInfo();
1127  G4cout << "Position:" << G4endl << G4endl ;
1128  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1129  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1130  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1131  G4cout.precision(oldprc) ;
1132  G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
1133  JustWarning, "Point p is outside !?" );
1134  }
1135 #endif
1136 
1137  // Z planes
1138  //
1139  distz1=fDz-p.z();
1140  distz2=fDz+p.z();
1141  if (distz1<distz2)
1142  {
1143  safe=distz1;
1144  }
1145  else
1146  {
1147  safe=distz2;
1148  }
1149 
1150  trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
1151 
1152  // Transformed x into `box' system
1153  //
1154  cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
1155  disty1=(fDy-trany)*cosy;
1156  disty2=(fDy+trany)*cosy;
1157 
1158  if (disty1<safe) safe=disty1;
1159  if (disty2<safe) safe=disty2;
1160 
1161  tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
1162  cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
1163  distx1=(fDx-tranx)*cosx;
1164  distx2=(fDx+tranx)*cosx;
1165 
1166  if (distx1<safe) safe=distx1;
1167  if (distx2<safe) safe=distx2;
1168 
1169  if (safe<0) safe=0;
1170  return safe;
1171 }
1172 
1174 //
1175 // Create a List containing the transformed vertices
1176 // Ordering [0-3] -fDz cross section
1177 // [4-7] +fDz cross section such that [0] is below [4],
1178 // [1] below [5] etc.
1179 // Note:
1180 // Caller has deletion resposibility
1181 
1182 G4ThreeVectorList*
1183 G4Para::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1184 {
1185  G4ThreeVectorList *vertices;
1186  vertices=new G4ThreeVectorList();
1187  if (vertices)
1188  {
1189  vertices->reserve(8);
1190  G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
1191  -fDz*fTthetaSphi-fDy, -fDz);
1192  G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
1193  -fDz*fTthetaSphi-fDy, -fDz);
1194  G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
1195  -fDz*fTthetaSphi+fDy, -fDz);
1196  G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
1197  -fDz*fTthetaSphi+fDy, -fDz);
1198  G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
1199  +fDz*fTthetaSphi-fDy, +fDz);
1200  G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
1201  +fDz*fTthetaSphi-fDy, +fDz);
1202  G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
1203  +fDz*fTthetaSphi+fDy, +fDz);
1204  G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
1205  +fDz*fTthetaSphi+fDy, +fDz);
1206 
1207  vertices->push_back(pTransform.TransformPoint(vertex0));
1208  vertices->push_back(pTransform.TransformPoint(vertex1));
1209  vertices->push_back(pTransform.TransformPoint(vertex2));
1210  vertices->push_back(pTransform.TransformPoint(vertex3));
1211  vertices->push_back(pTransform.TransformPoint(vertex4));
1212  vertices->push_back(pTransform.TransformPoint(vertex5));
1213  vertices->push_back(pTransform.TransformPoint(vertex6));
1214  vertices->push_back(pTransform.TransformPoint(vertex7));
1215  }
1216  else
1217  {
1218  DumpInfo();
1219  G4Exception("G4Para::CreateRotatedVertices()",
1220  "GeomSolids0003", FatalException,
1221  "Error in allocation of vertices. Out of memory !");
1222  }
1223  return vertices;
1224 }
1225 
1227 //
1228 // GetEntityType
1229 
1230 G4GeometryType G4Para::GetEntityType() const
1231 {
1232  return G4String("G4Para");
1233 }
1234 
1236 //
1237 // Make a clone of the object
1238 //
1239 G4VSolid* G4Para::Clone() const
1240 {
1241  return new G4Para(*this);
1242 }
1243 
1245 //
1246 // Stream object contents to an output stream
1247 
1248 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
1249 {
1250  G4int oldprc = os.precision(16);
1251  os << "-----------------------------------------------------------\n"
1252  << " *** Dump for solid - " << GetName() << " ***\n"
1253  << " ===================================================\n"
1254  << " Solid type: G4Para\n"
1255  << " Parameters: \n"
1256  << " half length X: " << fDx/mm << " mm \n"
1257  << " half length Y: " << fDy/mm << " mm \n"
1258  << " half length Z: " << fDz/mm << " mm \n"
1259  << " std::tan(alpha) : " << fTalpha/degree << " degrees \n"
1260  << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree
1261  << " degrees \n"
1262  << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree
1263  << " degrees \n"
1264  << "-----------------------------------------------------------\n";
1265  os.precision(oldprc);
1266 
1267  return os;
1268 }
1269 
1271 //
1272 // GetPointOnPlane
1273 // Auxiliary method for Get Point on Surface
1274 //
1275 
1276 G4ThreeVector G4Para::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1277  G4ThreeVector p2, G4ThreeVector p3,
1278  G4double& area) const
1279 {
1280  G4double lambda1, lambda2, chose, aOne, aTwo;
1281  G4ThreeVector t, u, v, w, Area, normal;
1282 
1283  t = p1 - p0;
1284  u = p2 - p1;
1285  v = p3 - p2;
1286  w = p0 - p3;
1287 
1288  Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
1289  w.z()*v.x() - w.x()*v.z(),
1290  w.x()*v.y() - w.y()*v.x());
1291 
1292  aOne = 0.5*Area.mag();
1293 
1294  Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
1295  t.z()*u.x() - t.x()*u.z(),
1296  t.x()*u.y() - t.y()*u.x());
1297 
1298  aTwo = 0.5*Area.mag();
1299 
1300  area = aOne + aTwo;
1301 
1302  chose = RandFlat::shoot(0.,aOne+aTwo);
1303 
1304  if( (chose>=0.) && (chose < aOne) )
1305  {
1306  lambda1 = RandFlat::shoot(0.,1.);
1307  lambda2 = RandFlat::shoot(0.,lambda1);
1308  return (p2+lambda1*v+lambda2*w);
1309  }
1310 
1311  // else
1312 
1313  lambda1 = RandFlat::shoot(0.,1.);
1314  lambda2 = RandFlat::shoot(0.,lambda1);
1315  return (p0+lambda1*t+lambda2*u);
1316 }
1317 
1319 //
1320 // GetPointOnSurface
1321 //
1322 // Return a point (G4ThreeVector) randomly and uniformly
1323 // selected on the solid surface
1324 
1325 G4ThreeVector G4Para::GetPointOnSurface() const
1326 {
1327  G4ThreeVector One, Two, Three, Four, Five, Six;
1328  G4ThreeVector pt[8] ;
1329  G4double chose, aOne, aTwo, aThree, aFour, aFive, aSix;
1330 
1331  pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
1332  -fDz*fTthetaSphi-fDy, -fDz);
1333  pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
1334  -fDz*fTthetaSphi-fDy, -fDz);
1335  pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
1336  -fDz*fTthetaSphi+fDy, -fDz);
1337  pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
1338  -fDz*fTthetaSphi+fDy, -fDz);
1339  pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
1340  +fDz*fTthetaSphi-fDy, +fDz);
1341  pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
1342  +fDz*fTthetaSphi-fDy, +fDz);
1343  pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
1344  +fDz*fTthetaSphi+fDy, +fDz);
1345  pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
1346  +fDz*fTthetaSphi+fDy, +fDz);
1347 
1348  // make sure we provide the points in a clockwise fashion
1349 
1350  One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1351  Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1352  Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1353  Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1354  Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1355  Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1356 
1357  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
1358 
1359  if( (chose>=0.) && (chose<aOne) )
1360  { return One; }
1361  else if(chose>=aOne && chose<aOne+aTwo)
1362  { return Two; }
1363  else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
1364  { return Three; }
1365  else if(chose>=aOne+aTwo+aThree && chose<aOne+aTwo+aThree+aFour)
1366  { return Four; }
1367  else if(chose>=aOne+aTwo+aThree+aFour && chose<aOne+aTwo+aThree+aFour+aFive)
1368  { return Five; }
1369  return Six;
1370 }
1371 
1373 //
1374 // Methods for visualisation
1375 
1376 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1377 {
1378  scene.AddSolid (*this);
1379 }
1380 
1381 G4Polyhedron* G4Para::CreatePolyhedron () const
1382 {
1383  G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1384  G4double alpha = std::atan(fTalpha);
1385  G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1386  +fTthetaSphi*fTthetaSphi));
1387 
1388  return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
1389 }
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:869
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:67
G4Para::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Para.cc:1325
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:551
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:646
G4AffineTransform::IsRotated
G4bool IsRotated() const
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:1248
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:1381
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:1230
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:1376
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:110
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
FatalException
Definition: G4ExceptionSeverity.hh:60
G4Para::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Para.cc:1183
degree
static const double degree
Definition: G4SIunits.hh:125
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4Para::Clone
G4VSolid * Clone() const
Definition: G4Para.cc:1239
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:102
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:1276