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