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