source: trunk/source/geometry/solids/CSG/src/G4Para.cc@ 1349

Last change on this file since 1349 was 1337, checked in by garnier, 15 years ago

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