source: trunk/source/geometry/solids/specific/src/G4VTwistedFaceted.cc@ 1127

Last change on this file since 1127 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 39.1 KB
RevLine 
[831]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// $Id: G4VTwistedFaceted.cc,v 1.18 2007/05/25 09:42:34 gcosmo Exp $
[1058]27// GEANT4 tag $Name: geant4-09-02-ref-02 $
[831]28//
29//
30// --------------------------------------------------------------------
31// GEANT 4 class source file
32//
33//
34// G4VTwistedFaceted.cc
35//
36// Author:
37//
38// 04-Nov-2004 - O.Link (Oliver.Link@cern.ch)
39//
40// --------------------------------------------------------------------
41
42#include "G4VTwistedFaceted.hh"
43
44#include "G4VoxelLimits.hh"
45#include "G4AffineTransform.hh"
46#include "G4SolidExtentList.hh"
47#include "G4ClippablePolygon.hh"
48#include "G4VPVParameterisation.hh"
49#include "G4GeometryTolerance.hh"
50#include "meshdefs.hh"
51
52#include "G4VGraphicsScene.hh"
53#include "G4Polyhedron.hh"
54#include "G4VisExtent.hh"
55#include "G4NURBS.hh"
56#include "G4NURBStube.hh"
57#include "G4NURBScylinder.hh"
58#include "G4NURBStubesector.hh"
59
60#include "Randomize.hh"
61
62//=====================================================================
63//* constructors ------------------------------------------------------
64
65G4VTwistedFaceted::
66G4VTwistedFaceted( const G4String &pname, // Name of instance
67 G4double PhiTwist, // twist angle
68 G4double pDz, // half z length
69 G4double pTheta, // direction between end planes
70 G4double pPhi, // defined by polar and azim. angles
71 G4double pDy1, // half y length at -pDz
72 G4double pDx1, // half x length at -pDz,-pDy
73 G4double pDx2, // half x length at -pDz,+pDy
74 G4double pDy2, // half y length at +pDz
75 G4double pDx3, // half x length at +pDz,-pDy
76 G4double pDx4, // half x length at +pDz,+pDy
77 G4double pAlph // tilt angle
78 )
79 : G4VSolid(pname),
80 fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
81 fSide90(0), fSide180(0), fSide270(0),
82 fSurfaceArea(0.), fpPolyhedron(0)
83{
84
85 G4double pDytmp ;
86 G4double fDxUp ;
87 G4double fDxDown ;
88
89 fDx1 = pDx1 ;
90 fDx2 = pDx2 ;
91 fDx3 = pDx3 ;
92 fDx4 = pDx4 ;
93 fDy1 = pDy1 ;
94 fDy2 = pDy2 ;
95 fDz = pDz ;
96
97 G4double kAngTolerance
98 = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
99
100 // maximum values
101 //
102 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
103 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
104 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
105 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
106
107 // planarity check
108 //
109 if ( fDx1 != fDx2 && fDx3 != fDx4 )
110 {
111 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
112 if ( std::fabs(pDytmp - fDy2) > kCarTolerance )
113 {
114 G4cerr << "ERROR - G4VTwistedFaceted::G4VTwistedFaceted(): "
115 << GetName() << G4endl
116 << " Not planar ! - " << G4endl
117 << "fDy2 is " << fDy2 << " but should be "
118 << pDytmp << "." << G4endl ;
119 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "InvalidSetup",
120 FatalException, "Not planar surface in untwisted Trapezoid.");
121 }
122 }
123
124#ifdef G4TWISTDEBUG
125 if ( fDx1 == fDx2 && fDx3 == fDx4 )
126 {
127 G4cout << "Trapezoid is a box" << G4endl ;
128 }
129
130#endif
131
132 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
133 {
134 G4cerr << "ERROR - G4VTwistedFaceted::G4VTwistedFaceted(): "
135 << GetName() << G4endl
136 << " Not planar ! - " << G4endl
137 << "One endcap is rectengular, the other is a trapezoid." << G4endl
138 << "For planarity reasons they have to be rectangles or trapezoids "
139 << G4endl
140 << "on both sides."
141 << G4endl ;
142 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "InvalidSetup",
143 FatalException, "Not planar surface in untwisted Trapezoid.");
144 }
145
146 // twist angle
147 //
148 fPhiTwist = PhiTwist ;
149
150 // tilt angle
151 //
152 fAlph = pAlph ;
153 fTAlph = std::tan(fAlph) ;
154
155 fTheta = pTheta ;
156 fPhi = pPhi ;
157
158 // dx in surface equation
159 //
160 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
161
162 // dy in surface equation
163 //
164 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
165
166 if ( ! ( ( fDx1 > 2*kCarTolerance)
167 && ( fDx2 > 2*kCarTolerance)
168 && ( fDx3 > 2*kCarTolerance)
169 && ( fDx4 > 2*kCarTolerance)
170 && ( fDy1 > 2*kCarTolerance)
171 && ( fDy2 > 2*kCarTolerance)
172 && ( fDz > 2*kCarTolerance)
173 && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
174 && ( std::fabs(fPhiTwist) < pi/2 )
175 && ( std::fabs(fAlph) < pi/2 )
176 && ( fTheta < pi/2 && fTheta >= 0 ) )
177 )
178 {
179 G4cerr << "ERROR - G4VTwistedFaceted()::G4VTwistedFaceted(): "
180 << GetName() << G4endl
181 << " Dimensions too small or too big! - " << G4endl
182 << "fDx 1-4 = " << fDx1/cm << ", " << fDx2/cm << ", "
183 << fDx3/cm << ", " << fDx4/cm << " cm" << G4endl
184 << "fDy 1-2 = " << fDy1/cm << ", " << fDy2/cm << ", "
185 << " cm" << G4endl
186 << "fDz = " << fDz/cm << " cm" << G4endl
187 << " twistangle " << fPhiTwist/deg << " deg" << G4endl
188 << " phi,theta = " << fPhi/deg << ", " << fTheta/deg
189 << " deg" << G4endl ;
190 G4Exception("G4TwistedTrap::G4VTwistedFaceted()",
191 "InvalidSetup", FatalException,
192 "Invalid dimensions. Too small, or twist angle too big.");
193 }
194 CreateSurfaces();
195 fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
196}
197
198
199//=====================================================================
200//* Fake default constructor ------------------------------------------
201
202G4VTwistedFaceted::G4VTwistedFaceted( __void__& a )
203 : G4VSolid(a),
204 fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
205 fSide90(0), fSide180(0), fSide270(0),
206 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
207{
208}
209
210//=====================================================================
211//* destructor --------------------------------------------------------
212
213G4VTwistedFaceted::~G4VTwistedFaceted()
214{
215 if (fLowerEndcap) delete fLowerEndcap ;
216 if (fUpperEndcap) delete fUpperEndcap ;
217
218 if (fSide0) delete fSide0 ;
219 if (fSide90) delete fSide90 ;
220 if (fSide180) delete fSide180 ;
221 if (fSide270) delete fSide270 ;
222 if (fpPolyhedron) delete fpPolyhedron;
223}
224
225//=====================================================================
226//* ComputeDimensions -------------------------------------------------
227
228void G4VTwistedFaceted::ComputeDimensions(G4VPVParameterisation* ,
229 const G4int ,
230 const G4VPhysicalVolume* )
231{
232 G4Exception("G4VTwistedFaceted::ComputeDimensions()",
233 "NotSupported", FatalException,
234 "G4VTwistedFaceted does not support Parameterisation.");
235}
236
237//=====================================================================
238//* CalculateExtent ---------------------------------------------------
239
240G4bool
241G4VTwistedFaceted::CalculateExtent( const EAxis pAxis,
242 const G4VoxelLimits &pVoxelLimit,
243 const G4AffineTransform &pTransform,
244 G4double &pMin,
245 G4double &pMax ) const
246{
247 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
248
249 if (!pTransform.IsRotated())
250 {
251 // Special case handling for unrotated boxes
252 // Compute x/y/z mins and maxs respecting limits, with early returns
253 // if outside limits. Then switch() on pAxis
254
255 G4double xoffset,xMin,xMax;
256 G4double yoffset,yMin,yMax;
257 G4double zoffset,zMin,zMax;
258
259 xoffset = pTransform.NetTranslation().x() ;
260 xMin = xoffset - maxRad ;
261 xMax = xoffset + maxRad ;
262
263 if (pVoxelLimit.IsXLimited())
264 {
265 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
266 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false;
267 else
268 {
269 if (xMin < pVoxelLimit.GetMinXExtent())
270 {
271 xMin = pVoxelLimit.GetMinXExtent() ;
272 }
273 if (xMax > pVoxelLimit.GetMaxXExtent())
274 {
275 xMax = pVoxelLimit.GetMaxXExtent() ;
276 }
277 }
278 }
279 yoffset = pTransform.NetTranslation().y() ;
280 yMin = yoffset - maxRad ;
281 yMax = yoffset + maxRad ;
282
283 if (pVoxelLimit.IsYLimited())
284 {
285 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
286 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false;
287 else
288 {
289 if (yMin < pVoxelLimit.GetMinYExtent())
290 {
291 yMin = pVoxelLimit.GetMinYExtent() ;
292 }
293 if (yMax > pVoxelLimit.GetMaxYExtent())
294 {
295 yMax = pVoxelLimit.GetMaxYExtent() ;
296 }
297 }
298 }
299 zoffset = pTransform.NetTranslation().z() ;
300 zMin = zoffset - fDz ;
301 zMax = zoffset + fDz ;
302
303 if (pVoxelLimit.IsZLimited())
304 {
305 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
306 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false;
307 else
308 {
309 if (zMin < pVoxelLimit.GetMinZExtent())
310 {
311 zMin = pVoxelLimit.GetMinZExtent() ;
312 }
313 if (zMax > pVoxelLimit.GetMaxZExtent())
314 {
315 zMax = pVoxelLimit.GetMaxZExtent() ;
316 }
317 }
318 }
319 switch (pAxis)
320 {
321 case kXAxis:
322 pMin = xMin ;
323 pMax = xMax ;
324 break ;
325 case kYAxis:
326 pMin=yMin;
327 pMax=yMax;
328 break;
329 case kZAxis:
330 pMin=zMin;
331 pMax=zMax;
332 break;
333 default:
334 break;
335 }
336 pMin -= kCarTolerance ;
337 pMax += kCarTolerance ;
338
339 return true;
340 }
341 else // General rotated case - create and clip mesh to boundaries
342 {
343 G4bool existsAfterClip = false ;
344 G4ThreeVectorList* vertices ;
345
346 pMin = +kInfinity ;
347 pMax = -kInfinity ;
348
349 // Calculate rotated vertex coordinates
350
351 vertices = CreateRotatedVertices(pTransform) ;
352 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
353 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
354 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
355
356 if (pVoxelLimit.IsLimited(pAxis) == false)
357 {
358 if ( pMin != kInfinity || pMax != -kInfinity )
359 {
360 existsAfterClip = true ;
361
362 // Add 2*tolerance to avoid precision troubles
363
364 pMin -= kCarTolerance;
365 pMax += kCarTolerance;
366 }
367 }
368 else
369 {
370 G4ThreeVector clipCentre(
371 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
372 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
373 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
374
375 if ( pMin != kInfinity || pMax != -kInfinity )
376 {
377 existsAfterClip = true ;
378
379 // Check to see if endpoints are in the solid
380
381 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
382
383 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
384 != kOutside)
385 {
386 pMin = pVoxelLimit.GetMinExtent(pAxis);
387 }
388 else
389 {
390 pMin -= kCarTolerance;
391 }
392 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
393
394 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
395 != kOutside)
396 {
397 pMax = pVoxelLimit.GetMaxExtent(pAxis);
398 }
399 else
400 {
401 pMax += kCarTolerance;
402 }
403 }
404
405 // Check for case where completely enveloping clipping volume
406 // If point inside then we are confident that the solid completely
407 // envelopes the clipping volume. Hence set min/max extents according
408 // to clipping volume extents along the specified axis.
409
410 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
411 != kOutside)
412 {
413 existsAfterClip = true ;
414 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
415 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
416 }
417 }
418 delete vertices;
419 return existsAfterClip;
420 }
421
422
423}
424
425G4ThreeVectorList* G4VTwistedFaceted::
426CreateRotatedVertices(const G4AffineTransform& pTransform) const
427{
428
429 G4ThreeVectorList* vertices = new G4ThreeVectorList();
430 vertices->reserve(8);
431
432 if (vertices)
433 {
434
435 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
436
437 G4ThreeVector vertex0(-maxRad,-maxRad,-fDz) ;
438 G4ThreeVector vertex1(maxRad,-maxRad,-fDz) ;
439 G4ThreeVector vertex2(maxRad,maxRad,-fDz) ;
440 G4ThreeVector vertex3(-maxRad,maxRad,-fDz) ;
441 G4ThreeVector vertex4(-maxRad,-maxRad,fDz) ;
442 G4ThreeVector vertex5(maxRad,-maxRad,fDz) ;
443 G4ThreeVector vertex6(maxRad,maxRad,fDz) ;
444 G4ThreeVector vertex7(-maxRad,maxRad,fDz) ;
445
446 vertices->push_back(pTransform.TransformPoint(vertex0));
447 vertices->push_back(pTransform.TransformPoint(vertex1));
448 vertices->push_back(pTransform.TransformPoint(vertex2));
449 vertices->push_back(pTransform.TransformPoint(vertex3));
450 vertices->push_back(pTransform.TransformPoint(vertex4));
451 vertices->push_back(pTransform.TransformPoint(vertex5));
452 vertices->push_back(pTransform.TransformPoint(vertex6));
453 vertices->push_back(pTransform.TransformPoint(vertex7));
454 }
455 else
456 {
457 DumpInfo();
458 G4Exception("G4VTwistedFaceted::CreateRotatedVertices()",
459 "FatalError", FatalException,
460 "Error in allocation of vertices. Out of memory !");
461 }
462 return vertices;
463}
464
465//=====================================================================
466//* Inside ------------------------------------------------------------
467
468EInside G4VTwistedFaceted::Inside(const G4ThreeVector& p) const
469{
470
471 G4ThreeVector *tmpp;
472 EInside *tmpin;
473 if (fLastInside.p == p) {
474 return fLastInside.inside;
475 } else {
476 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
477 tmpin = const_cast<EInside*>(&(fLastInside.inside));
478 tmpp->set(p.x(), p.y(), p.z());
479 }
480
481 *tmpin = kOutside ;
482
483 G4double phi = p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0
484 G4double cphi = std::cos(-phi) ;
485 G4double sphi = std::sin(-phi) ;
486
487 G4double px = p.x() + fdeltaX * ( -phi/fPhiTwist) ; // shift
488 G4double py = p.y() + fdeltaY * ( -phi/fPhiTwist) ;
489 G4double pz = p.z() ;
490
491 G4double posx = px * cphi - py * sphi ; // rotation
492 G4double posy = px * sphi + py * cphi ;
493 G4double posz = pz ;
494
495 G4double xMin = Xcoef(posy,phi,fTAlph) - 2*Xcoef(posy,phi,0.) ;
496 G4double xMax = Xcoef(posy,phi,fTAlph) ;
497
498 G4double yMax = GetValueB(phi)/2. ; // b(phi)/2 is limit
499 G4double yMin = -yMax ;
500
501#ifdef G4TWISTDEBUG
502
503 G4cout << "inside called: p = " << p << G4endl ;
504 G4cout << "fDx1 = " << fDx1 << G4endl ;
505 G4cout << "fDx2 = " << fDx2 << G4endl ;
506 G4cout << "fDx3 = " << fDx3 << G4endl ;
507 G4cout << "fDx4 = " << fDx4 << G4endl ;
508
509 G4cout << "fDy1 = " << fDy1 << G4endl ;
510 G4cout << "fDy2 = " << fDy2 << G4endl ;
511
512 G4cout << "fDz = " << fDz << G4endl ;
513
514 G4cout << "Tilt angle alpha = " << fAlph << G4endl ;
515 G4cout << "phi,theta = " << fPhi << " , " << fTheta << G4endl ;
516
517 G4cout << "Twist angle = " << fPhiTwist << G4endl ;
518
519 G4cout << "posx = " << posx << G4endl ;
520 G4cout << "posy = " << posy << G4endl ;
521 G4cout << "xMin = " << xMin << G4endl ;
522 G4cout << "xMax = " << xMax << G4endl ;
523 G4cout << "yMin = " << yMin << G4endl ;
524 G4cout << "yMax = " << yMax << G4endl ;
525
526#endif
527
528
529 if ( posx <= xMax - kCarTolerance*0.5
530 && posx >= xMin + kCarTolerance*0.5 )
531 {
532 if ( posy <= yMax - kCarTolerance*0.5
533 && posy >= yMin + kCarTolerance*0.5 )
534 {
535 if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) *tmpin = kInside ;
536 else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
537 }
538 else if ( posy <= yMax + kCarTolerance*0.5
539 && posy >= yMin - kCarTolerance*0.5 )
540 {
541 if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
542 }
543 }
544 else if ( posx <= xMax + kCarTolerance*0.5
545 && posx >= xMin - kCarTolerance*0.5 )
546 {
547 if ( posy <= yMax + kCarTolerance*0.5
548 && posy >= yMin - kCarTolerance*0.5 )
549 {
550 if (std::fabs(posz) <= fDz + kCarTolerance*0.5) *tmpin = kSurface ;
551 }
552 }
553
554#ifdef G4TWISTDEBUG
555 G4cout << "inside = " << fLastInside.inside << G4endl ;
556#endif
557
558 return fLastInside.inside;
559
560}
561
562//=====================================================================
563//* SurfaceNormal -----------------------------------------------------
564
565G4ThreeVector G4VTwistedFaceted::SurfaceNormal(const G4ThreeVector& p) const
566{
567 //
568 // return the normal unit vector to the Hyperbolical Surface at a point
569 // p on (or nearly on) the surface
570 //
571 // Which of the three or four surfaces are we closest to?
572 //
573
574 if (fLastNormal.p == p)
575 {
576 return fLastNormal.vec;
577 }
578
579 G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
580 G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
581 G4VTwistSurface **tmpsurface = const_cast<G4VTwistSurface**>(fLastNormal.surface);
582 tmpp->set(p.x(), p.y(), p.z());
583
584 G4double distance = kInfinity;
585
586 G4VTwistSurface *surfaces[6];
587
588 surfaces[0] = fSide0 ;
589 surfaces[1] = fSide90 ;
590 surfaces[2] = fSide180 ;
591 surfaces[3] = fSide270 ;
592 surfaces[4] = fLowerEndcap;
593 surfaces[5] = fUpperEndcap;
594
595 G4ThreeVector xx;
596 G4ThreeVector bestxx;
597 G4int i;
598 G4int besti = -1;
599 for (i=0; i< 6; i++)
600 {
601 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
602 if (tmpdistance < distance)
603 {
604 distance = tmpdistance;
605 bestxx = xx;
606 besti = i;
607 }
608 }
609
610 tmpsurface[0] = surfaces[besti];
611 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
612
613 return fLastNormal.vec;
614}
615
616//=====================================================================
617//* DistanceToIn (p, v) -----------------------------------------------
618
619G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p,
620 const G4ThreeVector& v ) const
621{
622
623 // DistanceToIn (p, v):
624 // Calculate distance to surface of shape from `outside'
625 // along with the v, allowing for tolerance.
626 // The function returns kInfinity if no intersection or
627 // just grazing within tolerance.
628
629 //
630 // checking last value
631 //
632
633 G4ThreeVector *tmpp;
634 G4ThreeVector *tmpv;
635 G4double *tmpdist;
636 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
637 {
638 return fLastDistanceToIn.value;
639 }
640 else
641 {
642 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
643 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
644 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
645 tmpp->set(p.x(), p.y(), p.z());
646 tmpv->set(v.x(), v.y(), v.z());
647 }
648
649 //
650 // Calculate DistanceToIn(p,v)
651 //
652
653 EInside currentside = Inside(p);
654
655 if (currentside == kInside)
656 {
657 }
658 else if (currentside == kSurface)
659 {
660 // particle is just on a boundary.
661 // if the particle is entering to the volume, return 0
662 //
663 G4ThreeVector normal = SurfaceNormal(p);
664 if (normal*v < 0)
665 {
666 *tmpdist = 0;
667 return fLastDistanceToInWithV.value;
668 }
669 }
670
671 // now, we can take smallest positive distance.
672
673 // Initialize
674 //
675 G4double distance = kInfinity;
676
677 // Find intersections and choose nearest one
678 //
679 G4VTwistSurface *surfaces[6];
680
681 surfaces[0] = fSide0;
682 surfaces[1] = fSide90 ;
683 surfaces[2] = fSide180 ;
684 surfaces[3] = fSide270 ;
685 surfaces[4] = fLowerEndcap;
686 surfaces[5] = fUpperEndcap;
687
688 G4ThreeVector xx;
689 G4ThreeVector bestxx;
690 G4int i;
691 G4int besti = -1;
692 for (i=0; i < 6 ; i++)
693 //for (i=1; i < 2 ; i++)
694 {
695
696#ifdef G4TWISTDEBUG
697 G4cout << G4endl << "surface " << i << ": " << G4endl << G4endl ;
698#endif
699 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
700#ifdef G4TWISTDEBUG
701 G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ;
702 G4cout << "intersection point = " << xx << G4endl ;
703#endif
704 if (tmpdistance < distance)
705 {
706 distance = tmpdistance;
707 bestxx = xx;
708 besti = i;
709 }
710 }
711
712#ifdef G4TWISTDEBUG
713 G4cout << "best distance = " << distance << G4endl ;
714#endif
715
716 *tmpdist = distance;
717 // timer.Stop();
718 return fLastDistanceToInWithV.value;
719}
720
721
722G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p) const
723{
724 // DistanceToIn(p):
725 // Calculate distance to surface of shape from `outside',
726 // allowing for tolerance
727 //
728
729 //
730 // checking last value
731 //
732
733 G4ThreeVector *tmpp;
734 G4double *tmpdist;
735 if (fLastDistanceToIn.p == p)
736 {
737 return fLastDistanceToIn.value;
738 }
739 else
740 {
741 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
742 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
743 tmpp->set(p.x(), p.y(), p.z());
744 }
745
746 //
747 // Calculate DistanceToIn(p)
748 //
749
750 EInside currentside = Inside(p);
751
752 switch (currentside)
753 {
754 case (kInside) :
755 {
756 }
757
758 case (kSurface) :
759 {
760 *tmpdist = 0.;
761 return fLastDistanceToIn.value;
762 }
763
764 case (kOutside) :
765 {
766 // Initialize
767 //
768 G4double distance = kInfinity;
769
770 // Find intersections and choose nearest one
771 //
772 G4VTwistSurface *surfaces[6];
773
774 surfaces[0] = fSide0;
775 surfaces[1] = fSide90 ;
776 surfaces[2] = fSide180 ;
777 surfaces[3] = fSide270 ;
778 surfaces[4] = fLowerEndcap;
779 surfaces[5] = fUpperEndcap;
780
781 G4int i;
782 G4int besti = -1;
783 G4ThreeVector xx;
784 G4ThreeVector bestxx;
785 for (i=0; i< 6; i++)
786 {
787 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
788 if (tmpdistance < distance)
789 {
790 distance = tmpdistance;
791 bestxx = xx;
792 besti = i;
793 }
794 }
795 *tmpdist = distance;
796 return fLastDistanceToIn.value;
797 }
798
799 default :
800 {
801 G4Exception("G4VTwistedFaceted::DistanceToIn(p)", "InvalidCondition",
802 FatalException, "Unknown point location!");
803 }
804 } // switch end
805
806 return 0;
807}
808
809
810//=====================================================================
811//* DistanceToOut (p, v) ----------------------------------------------
812
813G4double
814G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p,
815 const G4ThreeVector& v,
816 const G4bool calcNorm,
817 G4bool *validNorm,
818 G4ThreeVector *norm ) const
819{
820 // DistanceToOut (p, v):
821 // Calculate distance to surface of shape from `inside'
822 // along with the v, allowing for tolerance.
823 // The function returns kInfinity if no intersection or
824 // just grazing within tolerance.
825
826 //
827 // checking last value
828 //
829
830 G4ThreeVector *tmpp;
831 G4ThreeVector *tmpv;
832 G4double *tmpdist;
833 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
834 {
835 return fLastDistanceToOutWithV.value;
836 }
837 else
838 {
839 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
840 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
841 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
842 tmpp->set(p.x(), p.y(), p.z());
843 tmpv->set(v.x(), v.y(), v.z());
844 }
845
846 //
847 // Calculate DistanceToOut(p,v)
848 //
849
850 EInside currentside = Inside(p);
851
852 if (currentside == kOutside)
853 {
854 }
855 else if (currentside == kSurface)
856 {
857 // particle is just on a boundary.
858 // if the particle is exiting from the volume, return 0
859 //
860 G4ThreeVector normal = SurfaceNormal(p);
861 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
862 if (normal*v > 0)
863 {
864 if (calcNorm)
865 {
866 *norm = (blockedsurface->GetNormal(p, true));
867 *validNorm = blockedsurface->IsValidNorm();
868 }
869 *tmpdist = 0.;
870 // timer.Stop();
871 return fLastDistanceToOutWithV.value;
872 }
873 }
874
875 // now, we can take smallest positive distance.
876
877 // Initialize
878 G4double distance = kInfinity;
879
880 // find intersections and choose nearest one.
881 G4VTwistSurface *surfaces[6];
882
883 surfaces[0] = fSide0;
884 surfaces[1] = fSide90 ;
885 surfaces[2] = fSide180 ;
886 surfaces[3] = fSide270 ;
887 surfaces[4] = fLowerEndcap;
888 surfaces[5] = fUpperEndcap;
889
890 G4int i;
891 G4int besti = -1;
892 G4ThreeVector xx;
893 G4ThreeVector bestxx;
894 for (i=0; i< 6 ; i++) {
895 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
896 if (tmpdistance < distance)
897 {
898 distance = tmpdistance;
899 bestxx = xx;
900 besti = i;
901 }
902 }
903
904 if (calcNorm)
905 {
906 if (besti != -1)
907 {
908 *norm = (surfaces[besti]->GetNormal(p, true));
909 *validNorm = surfaces[besti]->IsValidNorm();
910 }
911 }
912
913 *tmpdist = distance;
914 // timer.Stop();
915 return fLastDistanceToOutWithV.value;
916}
917
918
919G4double G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p ) const
920{
921 // DistanceToOut(p):
922 // Calculate distance to surface of shape from `inside',
923 // allowing for tolerance
924
925 //
926 // checking last value
927 //
928
929
930 G4ThreeVector *tmpp;
931 G4double *tmpdist;
932 if (fLastDistanceToOut.p == p)
933 {
934 return fLastDistanceToOut.value;
935 }
936 else
937 {
938 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
939 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
940 tmpp->set(p.x(), p.y(), p.z());
941 }
942
943 //
944 // Calculate DistanceToOut(p)
945 //
946
947 EInside currentside = Inside(p);
948
949 switch (currentside)
950 {
951 case (kOutside) :
952 {
953#ifdef G4SPECSDEBUG
954 G4cout.precision(16) ;
955 G4cout << G4endl ;
956 DumpInfo();
957 G4cout << "Position:" << G4endl << G4endl ;
958 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
959 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
960 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
961 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "Notification",
962 JustWarning, "Point p is outside !?" );
963#endif
964 }
965 case (kSurface) :
966 {
967 *tmpdist = 0.;
968 return fLastDistanceToOut.value;
969 }
970
971 case (kInside) :
972 {
973 // Initialize
974 //
975 G4double distance = kInfinity;
976
977 // find intersections and choose nearest one
978 //
979 G4VTwistSurface *surfaces[6];
980
981 surfaces[0] = fSide0;
982 surfaces[1] = fSide90 ;
983 surfaces[2] = fSide180 ;
984 surfaces[3] = fSide270 ;
985 surfaces[4] = fLowerEndcap;
986 surfaces[5] = fUpperEndcap;
987
988 G4int i;
989 G4int besti = -1;
990 G4ThreeVector xx;
991 G4ThreeVector bestxx;
992 for (i=0; i< 6; i++)
993 {
994 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
995 if (tmpdistance < distance)
996 {
997 distance = tmpdistance;
998 bestxx = xx;
999 besti = i;
1000 }
1001 }
1002 *tmpdist = distance;
1003
1004 return fLastDistanceToOut.value;
1005 }
1006
1007 default :
1008 {
1009 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "InvalidCondition",
1010 FatalException, "Unknown point location!");
1011 }
1012 } // switch end
1013
1014 return kInfinity;
1015}
1016
1017
1018//=====================================================================
1019//* StreamInfo --------------------------------------------------------
1020
1021std::ostream& G4VTwistedFaceted::StreamInfo(std::ostream& os) const
1022{
1023 //
1024 // Stream object contents to an output stream
1025 //
1026 os << "-----------------------------------------------------------\n"
1027 << " *** Dump for solid - " << GetName() << " ***\n"
1028 << " ===================================================\n"
1029 << " Solid type: G4VTwistedFaceted\n"
1030 << " Parameters: \n"
1031 << " polar angle theta = " << fTheta/degree << " deg" << G4endl
1032 << " azimuthal angle phi = " << fPhi/degree << " deg" << G4endl
1033 << " tilt angle alpha = " << fAlph/degree << " deg" << G4endl
1034 << " TWIST angle = " << fPhiTwist/degree << " deg" << G4endl
1035 << " Half length along y (lower endcap) = " << fDy1/cm << " cm"
1036 << G4endl
1037 << " Half length along x (lower endcap, bottom) = " << fDx1/cm << " cm"
1038 << G4endl
1039 << " Half length along x (lower endcap, top) = " << fDx2/cm << " cm"
1040 << G4endl
1041 << " Half length along y (upper endcap) = " << fDy2/cm << " cm"
1042 << G4endl
1043 << " Half length along x (upper endcap, bottom) = " << fDx3/cm << " cm"
1044 << G4endl
1045 << " Half length along x (upper endcap, top) = " << fDx4/cm << " cm"
1046 << G4endl
1047 << "-----------------------------------------------------------\n";
1048
1049 return os;
1050}
1051
1052
1053//=====================================================================
1054//* DiscribeYourselfTo ------------------------------------------------
1055
1056void G4VTwistedFaceted::DescribeYourselfTo (G4VGraphicsScene& scene) const
1057{
1058 scene.AddSolid (*this);
1059}
1060
1061//=====================================================================
1062//* GetExtent ---------------------------------------------------------
1063
1064G4VisExtent G4VTwistedFaceted::GetExtent() const
1065{
1066 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
1067
1068 return G4VisExtent(-maxRad, maxRad ,
1069 -maxRad, maxRad ,
1070 -fDz, fDz );
1071}
1072
1073
1074//=====================================================================
1075//* CreateNUBS --------------------------------------------------------
1076
1077G4NURBS* G4VTwistedFaceted::CreateNURBS () const
1078{
1079 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
1080
1081 return new G4NURBStube(maxRad, maxRad, fDz);
1082 // Tube for now!!!
1083}
1084
1085
1086//=====================================================================
1087//* CreateSurfaces ----------------------------------------------------
1088
1089void G4VTwistedFaceted::CreateSurfaces()
1090{
1091
1092 // create 6 surfaces of TwistedTub.
1093
1094 if ( fDx1 == fDx2 && fDx3 == fDx4 ) // special case : Box
1095 {
1096 fSide0 = new G4TwistBoxSide("0deg", fPhiTwist, fDz, fTheta, fPhi,
1097 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*deg);
1098 fSide180 = new G4TwistBoxSide("180deg", fPhiTwist, fDz, fTheta, fPhi+pi,
1099 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*deg);
1100 }
1101 else // default general case
1102 {
1103 fSide0 = new G4TwistTrapAlphaSide("0deg" ,fPhiTwist, fDz, fTheta,
1104 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
1105 fSide180 = new G4TwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
1106 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
1107 }
1108
1109 // create parallel sides
1110 //
1111 fSide90 = new G4TwistTrapParallelSide("90deg", fPhiTwist, fDz, fTheta,
1112 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
1113 fSide270 = new G4TwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
1114 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
1115
1116 // create endcaps
1117 //
1118 fUpperEndcap = new G4TwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
1119 fDz, fAlph, fPhi, fTheta, 1 );
1120 fLowerEndcap = new G4TwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
1121 fDz, fAlph, fPhi, fTheta, -1 );
1122
1123 // Set neighbour surfaces
1124
1125 fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
1126 fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
1127 fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
1128 fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
1129 fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1130 fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1131
1132}
1133
1134
1135//=====================================================================
1136//* GetEntityType -----------------------------------------------------
1137
1138G4GeometryType G4VTwistedFaceted::GetEntityType() const
1139{
1140 return G4String("G4VTwistedFaceted");
1141}
1142
1143
1144//=====================================================================
1145//* GetPolyhedron -----------------------------------------------------
1146
1147G4Polyhedron* G4VTwistedFaceted::GetPolyhedron() const
1148{
1149 if (!fpPolyhedron ||
1150 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1151 fpPolyhedron->GetNumberOfRotationSteps())
1152 {
1153 delete fpPolyhedron;
1154 fpPolyhedron = CreatePolyhedron();
1155 }
1156
1157 return fpPolyhedron;
1158}
1159
1160
1161//=====================================================================
1162//* GetPointInSolid ---------------------------------------------------
1163
1164G4ThreeVector G4VTwistedFaceted::GetPointInSolid(G4double z) const
1165{
1166
1167
1168 // this routine is only used for a test
1169 // can be deleted ...
1170
1171 if ( z == fDz ) z -= 0.1*fDz ;
1172 if ( z == -fDz ) z += 0.1*fDz ;
1173
1174 G4double phi = z/(2*fDz)*fPhiTwist ;
1175
1176 return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1177}
1178
1179
1180//=====================================================================
1181//* GetPointOnSurface -------------------------------------------------
1182
1183G4ThreeVector G4VTwistedFaceted::GetPointOnSurface() const
1184{
1185
1186 G4double phi = CLHEP::RandFlat::shoot(-fPhiTwist/2.,fPhiTwist/2.);
1187 G4double u , umin, umax ; // variable for twisted surfaces
1188 G4double y ; // variable for flat surface (top and bottom)
1189
1190 // Compute the areas. Attention: Only correct for trapezoids
1191 // where the twisting is done along the z-axis. In the general case
1192 // the computed surface area is more difficult. However this simplification
1193 // does not affect the tracking through the solid.
1194
1195 G4double a1 = fSide0->GetSurfaceArea();
1196 G4double a2 = fSide90->GetSurfaceArea();
1197 G4double a3 = fSide180->GetSurfaceArea() ;
1198 G4double a4 = fSide270->GetSurfaceArea() ;
1199 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1200 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1201
1202#ifdef G4TWISTDEBUG
1203 G4cout << "Surface 0 deg = " << a1 << G4endl ;
1204 G4cout << "Surface 90 deg = " << a2 << G4endl ;
1205 G4cout << "Surface 180 deg = " << a3 << G4endl ;
1206 G4cout << "Surface 270 deg = " << a4 << G4endl ;
1207 G4cout << "Surface Lower = " << a5 << G4endl ;
1208 G4cout << "Surface Upper = " << a6 << G4endl ;
1209#endif
1210
1211 G4double chose = CLHEP::RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1212
1213 if(chose < a1)
1214 {
1215
1216 umin = fSide0->GetBoundaryMin(phi) ;
1217 umax = fSide0->GetBoundaryMax(phi) ;
1218 u = CLHEP::RandFlat::shoot(umin,umax) ;
1219
1220 return fSide0->SurfacePoint(phi, u, true) ; // point on 0deg surface
1221 }
1222
1223 else if( (chose >= a1) && (chose < a1 + a2 ) )
1224 {
1225
1226 umin = fSide90->GetBoundaryMin(phi) ;
1227 umax = fSide90->GetBoundaryMax(phi) ;
1228
1229 u = CLHEP::RandFlat::shoot(umin,umax) ;
1230
1231 return fSide90->SurfacePoint(phi, u, true); // point on 90deg surface
1232 }
1233
1234 else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1235 {
1236
1237 umin = fSide180->GetBoundaryMin(phi) ;
1238 umax = fSide180->GetBoundaryMax(phi) ;
1239 u = CLHEP::RandFlat::shoot(umin,umax) ;
1240
1241 return fSide180->SurfacePoint(phi, u, true); // point on 180 deg surface
1242 }
1243
1244 else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1245 {
1246
1247 umin = fSide270->GetBoundaryMin(phi) ;
1248 umax = fSide270->GetBoundaryMax(phi) ;
1249 u = CLHEP::RandFlat::shoot(umin,umax) ;
1250
1251 return fSide270->SurfacePoint(phi, u, true); // point on 270 deg surface
1252 }
1253
1254 else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1255 {
1256
1257 y = CLHEP::RandFlat::shoot(-fDy1,fDy1) ;
1258 umin = fLowerEndcap->GetBoundaryMin(y) ;
1259 umax = fLowerEndcap->GetBoundaryMax(y) ;
1260 u = CLHEP::RandFlat::shoot(umin,umax) ;
1261
1262 return fLowerEndcap->SurfacePoint(u,y,true); // point on lower endcap
1263 }
1264 else {
1265
1266 y = CLHEP::RandFlat::shoot(-fDy2,fDy2) ;
1267 umin = fUpperEndcap->GetBoundaryMin(y) ;
1268 umax = fUpperEndcap->GetBoundaryMax(y) ;
1269 u = CLHEP::RandFlat::shoot(umin,umax) ;
1270
1271 return fUpperEndcap->SurfacePoint(u,y,true) ; // point on upper endcap
1272
1273 }
1274}
1275
1276
1277//=====================================================================
1278//* CreatePolyhedron --------------------------------------------------
1279
1280G4Polyhedron* G4VTwistedFaceted::CreatePolyhedron () const
1281{
1282 // number of meshes
1283 const G4int m =
1284 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1285 const G4int n = m;
1286
1287 const G4int nnodes = 4*(m-1)*(n-2) + 2*m*m ;
1288 const G4int nfaces = 4*(m-1)*(n-1) + 2*(m-1)*(m-1) ;
1289
1290 G4Polyhedron *ph=new G4Polyhedron;
1291 typedef G4double G4double3[3];
1292 typedef G4int G4int4[4];
1293 G4double3* xyz = new G4double3[nnodes]; // number of nodes
1294 G4int4* faces = new G4int4[nfaces] ; // number of faces
1295
1296 fLowerEndcap->GetFacets(m,m,xyz,faces,0) ;
1297 fUpperEndcap->GetFacets(m,m,xyz,faces,1) ;
1298 fSide270->GetFacets(m,n,xyz,faces,2) ;
1299 fSide0->GetFacets(m,n,xyz,faces,3) ;
1300 fSide90->GetFacets(m,n,xyz,faces,4) ;
1301 fSide180->GetFacets(m,n,xyz,faces,5) ;
1302
1303 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1304
1305 return ph;
1306}
Note: See TracBrowser for help on using the repository browser.