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