Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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$
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#include "G4NURBS.hh"
57#include "G4NURBSbox.hh"
58
59using namespace CLHEP;
60
61// Private enum: Not for external use
62
63enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
64
65// used internally for normal routine
66
67enum ENSide {kNZ,kNX,kNY};
68
69/////////////////////////////////////////////////////////////////////
70//
71// Constructor - check and set half-widths
72
73void 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
99///////////////////////////////////////////////////////////////////////////
100//
101
102G4Para::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
119////////////////////////////////////////////////////////////////////////
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
125G4Para::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
154///////////////////////////////////////////////////////////////////////
155//
156// Fake default constructor - sets only member data and allocates memory
157// for usage restricted to object persistency.
158//
159G4Para::G4Para( __void__& a )
160 : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.),
161 fTalpha(0.), fTthetaCphi(0.), fTthetaSphi(0.)
162{
163}
164
165//////////////////////////////////////////////////////////////////////////
166//
167
168G4Para::~G4Para()
169{
170}
171
172//////////////////////////////////////////////////////////////////////////
173//
174// Copy constructor
175
176G4Para::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
183//////////////////////////////////////////////////////////////////////////
184//
185// Assignment operator
186
187G4Para& 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
206//////////////////////////////////////////////////////////////////////////
207//
208// Dispatch to parameterisation for replication mechanism dimension
209// computation & modification.
210
211void G4Para::ComputeDimensions( G4VPVParameterisation* p,
212 const G4int n,
213 const G4VPhysicalVolume* pRep )
214{
215 p->ComputeDimensions(*this,n,pRep);
216}
217
218
219//////////////////////////////////////////////////////////////
220//
221// Calculate extent under transform and specified limit
222
223G4bool 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
439/////////////////////////////////////////////////////////////////////////////
440//
441// Check in p is inside/on surface/outside solid
442
443EInside 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
477///////////////////////////////////////////////////////////////////////////
478//
479// Calculate side nearest to p, and return normal
480// If 2+ sides equidistant, first side's normal returned (arbitrarily)
481
482G4ThreeVector 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
551////////////////////////////////////////////////////////////////////////
552//
553// Algorithm for SurfaceNormal() following the original specification
554// for points not on the surface
555
556G4ThreeVector 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
635//////////////////////////////////////////////////////////////////////////////
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
651G4double 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
822////////////////////////////////////////////////////////////////////////////
823//
824// Calculate exact shortest distance to any boundary from outside
825// - Returns 0 is point inside
826
827G4double 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
869//////////////////////////////////////////////////////////////////////////
870//
871// Calculate distance to surface of shape from inside
872// Calculate distance to x/y/z planes - smallest is exiting distance
873
874G4double 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
1115/////////////////////////////////////////////////////////////////////////////
1116//
1117// Calculate exact shortest distance to any boundary from inside
1118// - Returns 0 is point outside
1119
1120G4double 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
1178////////////////////////////////////////////////////////////////////////////////
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
1187G4ThreeVectorList*
1188G4Para::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
1231//////////////////////////////////////////////////////////////////////////
1232//
1233// GetEntityType
1234
1235G4GeometryType G4Para::GetEntityType() const
1236{
1237 return G4String("G4Para");
1238}
1239
1240//////////////////////////////////////////////////////////////////////////
1241//
1242// Make a clone of the object
1243//
1244G4VSolid* G4Para::Clone() const
1245{
1246 return new G4Para(*this);
1247}
1248
1249//////////////////////////////////////////////////////////////////////////
1250//
1251// Stream object contents to an output stream
1252
1253std::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
1275//////////////////////////////////////////////////////////////////////////////
1276//
1277// GetPointOnPlane
1278// Auxiliary method for Get Point on Surface
1279//
1280
1281G4ThreeVector 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
1323/////////////////////////////////////////////////////////////////////////
1324//
1325// GetPointOnSurface
1326//
1327// Return a point (G4ThreeVector) randomly and uniformly
1328// selected on the solid surface
1329
1330G4ThreeVector 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
1377////////////////////////////////////////////////////////////////////////////
1378//
1379// Methods for visualisation
1380
1381void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1382{
1383 scene.AddSolid (*this);
1384}
1385
1386G4Polyhedron* 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
1396G4NURBS* G4Para::CreateNURBS () const
1397{
1398 // return new G4NURBSbox (fDx, fDy, fDz);
1399 return 0 ;
1400}
ESide
ESide
Definition: G4Cons.cc:68
kPZ
@ kPZ
Definition: G4Cons.cc:68
kMZ
@ kMZ
Definition: G4Cons.cc:68
kNZ
@ kNZ
Definition: G4Cons.cc:72
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
kPZ
@ kPZ
Definition: G4Para.cc:63
kPY
@ kPY
Definition: G4Para.cc:63
kMY
@ kMY
Definition: G4Para.cc:63
kPX
@ kPX
Definition: G4Para.cc:63
kUndef
@ kUndef
Definition: G4Para.cc:63
kMZ
@ kMZ
Definition: G4Para.cc:63
kMX
@ kMX
Definition: G4Para.cc:63
ENSide
ENSide
Definition: G4Para.cc:67
kNX
@ kNX
Definition: G4Para.cc:67
kNZ
@ kNZ
Definition: G4Para.cc:67
kNY
@ kNY
Definition: G4Para.cc:67
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:59
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4CSGSolid::fpPolyhedron
G4Polyhedron * fpPolyhedron
Definition: G4CSGSolid.hh:80
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
G4Para
Definition: G4Para.hh:77
G4Para::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Para.cc:1235
G4Para::Clone
G4VSolid * Clone() const
Definition: G4Para.cc:1244
G4Para::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Para.cc:1188
G4Para::CreateNURBS
G4NURBS * CreateNURBS() const
Definition: G4Para.cc:1396
G4Para::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Para.cc:1386
G4Para::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Para.cc:482
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:874
G4Para::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Para.cc:211
G4Para::operator=
G4Para & operator=(const G4Para &rhs)
Definition: G4Para.cc:187
G4Para::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4Para.cc:223
G4Para::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Para.cc:1253
G4Para::G4Para
G4Para(const G4String &pName, G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:102
G4Para::SetAllParameters
void SetAllParameters(G4double pDx, G4double pDy, G4double pDz, G4double pAlpha, G4double pTheta, G4double pPhi)
Definition: G4Para.cc:73
G4Para::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Para.cc:651
G4Para::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Para.cc:1381
G4Para::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Para.cc:443
G4Para::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Para.cc:1330
G4Para::~G4Para
virtual ~G4Para()
Definition: G4Para.cc:168
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:94
G4VSolid::GetName
G4String GetName() const
G4VSolid::ClipBetweenSections
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:376
G4VSolid::DumpInfo
void DumpInfo() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:345
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
EAxis
EAxis
Definition: geomdefs.hh:54
kYAxis
@ kYAxis
Definition: geomdefs.hh:54
kXAxis
@ kXAxis
Definition: geomdefs.hh:54
kZAxis
@ kZAxis
Definition: geomdefs.hh:54
EInside
EInside
Definition: geomdefs.hh:58
kInside
@ kInside
Definition: geomdefs.hh:58
kOutside
@ kOutside
Definition: geomdefs.hh:58
kSurface
@ kSurface
Definition: geomdefs.hh:58
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP
Definition: DoubConv.h:17