source: trunk/source/geometry/solids/specific/src/G4TwistedTubs.cc@ 1202

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

file release beta
File size: 39.3 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: G4TwistedTubs.cc,v 1.24 2007/05/18 07:39:56 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTubsSide.cc
36//
37// Author:
38// 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
39//
40// History:
41// 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
42// from original version in Jupiter-2.5.02 application.
43// --------------------------------------------------------------------
44
45#include "G4TwistedTubs.hh"
46
47#include "G4GeometryTolerance.hh"
48#include "G4VoxelLimits.hh"
49#include "G4AffineTransform.hh"
50#include "G4SolidExtentList.hh"
51#include "G4ClippablePolygon.hh"
52#include "G4VPVParameterisation.hh"
53#include "meshdefs.hh"
54
55#include "G4VGraphicsScene.hh"
56#include "G4Polyhedron.hh"
57#include "G4VisExtent.hh"
58#include "G4NURBS.hh"
59#include "G4NURBStube.hh"
60#include "G4NURBScylinder.hh"
61#include "G4NURBStubesector.hh"
62
63#include "Randomize.hh"
64
65//=====================================================================
66//* constructors ------------------------------------------------------
67
68G4TwistedTubs::G4TwistedTubs(const G4String &pname,
69 G4double twistedangle,
70 G4double endinnerrad,
71 G4double endouterrad,
72 G4double halfzlen,
73 G4double dphi)
74 : G4VSolid(pname), fDPhi(dphi),
75 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
76 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
77 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
78{
79 if (endinnerrad < DBL_MIN)
80 {
81 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
82 FatalException, "Invalid end-inner-radius!");
83 }
84
85 G4double sinhalftwist = std::sin(0.5 * twistedangle);
86
87 G4double endinnerradX = endinnerrad * sinhalftwist;
88 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
89 - endinnerradX * endinnerradX );
90
91 G4double endouterradX = endouterrad * sinhalftwist;
92 G4double outerrad = std::sqrt( endouterrad * endouterrad
93 - endouterradX * endouterradX );
94
95 // temporary treatment!!
96 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
97 CreateSurfaces();
98}
99
100G4TwistedTubs::G4TwistedTubs(const G4String &pname,
101 G4double twistedangle,
102 G4double endinnerrad,
103 G4double endouterrad,
104 G4double halfzlen,
105 G4int nseg,
106 G4double totphi)
107 : G4VSolid(pname),
108 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
109 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
110 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
111{
112
113 if (!nseg)
114 {
115 G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
116 << " Invalid nseg. nseg = " << nseg << G4endl;
117 }
118 if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
119 {
120 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
121 FatalException, "Invalid total-phi or end-inner-radius!");
122 }
123
124 G4double sinhalftwist = std::sin(0.5 * twistedangle);
125
126 G4double endinnerradX = endinnerrad * sinhalftwist;
127 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
128 - endinnerradX * endinnerradX );
129
130 G4double endouterradX = endouterrad * sinhalftwist;
131 G4double outerrad = std::sqrt( endouterrad * endouterrad
132 - endouterradX * endouterradX );
133
134 // temporary treatment!!
135 fDPhi = totphi / nseg;
136 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
137 CreateSurfaces();
138}
139
140G4TwistedTubs::G4TwistedTubs(const G4String &pname,
141 G4double twistedangle,
142 G4double innerrad,
143 G4double outerrad,
144 G4double negativeEndz,
145 G4double positiveEndz,
146 G4double dphi)
147 : G4VSolid(pname), fDPhi(dphi),
148 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
149 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
150 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
151{
152 if (innerrad < DBL_MIN)
153 {
154 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
155 FatalException, "Invalid end-inner-radius!");
156 }
157
158 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
159 CreateSurfaces();
160}
161
162G4TwistedTubs::G4TwistedTubs(const G4String &pname,
163 G4double twistedangle,
164 G4double innerrad,
165 G4double outerrad,
166 G4double negativeEndz,
167 G4double positiveEndz,
168 G4int nseg,
169 G4double totphi)
170 : G4VSolid(pname),
171 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
172 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
173 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
174{
175 if (!nseg)
176 {
177 G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
178 << " Invalid nseg. nseg = " << nseg << G4endl;
179 }
180 if (totphi == DBL_MIN || innerrad < DBL_MIN)
181 {
182 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
183 FatalException, "Invalid total-phi or end-inner-radius!");
184 }
185
186 fDPhi = totphi / nseg;
187 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
188 CreateSurfaces();
189}
190
191//=====================================================================
192//* Fake default constructor ------------------------------------------
193
194G4TwistedTubs::G4TwistedTubs( __void__& a )
195 : G4VSolid(a), fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
196 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), fCubicVolume(0.),
197 fSurfaceArea(0.), fpPolyhedron(0)
198{
199}
200
201//=====================================================================
202//* destructor --------------------------------------------------------
203
204G4TwistedTubs::~G4TwistedTubs()
205{
206 if (fLowerEndcap) { delete fLowerEndcap; }
207 if (fUpperEndcap) { delete fUpperEndcap; }
208 if (fLatterTwisted) { delete fLatterTwisted; }
209 if (fFormerTwisted) { delete fFormerTwisted; }
210 if (fInnerHype) { delete fInnerHype; }
211 if (fOuterHype) { delete fOuterHype; }
212 if (fpPolyhedron) { delete fpPolyhedron; }
213}
214
215//=====================================================================
216//* ComputeDimensions -------------------------------------------------
217
218void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
219 const G4int /* n */ ,
220 const G4VPhysicalVolume* /* pRep */ )
221{
222 G4Exception("G4TwistedTubs::ComputeDimensions()",
223 "NotSupported", FatalException,
224 "G4TwistedTubs does not support Parameterisation.");
225}
226
227
228//=====================================================================
229//* CalculateExtent ---------------------------------------------------
230
231G4bool G4TwistedTubs::CalculateExtent( const EAxis axis,
232 const G4VoxelLimits &voxelLimit,
233 const G4AffineTransform &transform,
234 G4double &min,
235 G4double &max ) const
236{
237
238 G4SolidExtentList extentList( axis, voxelLimit );
239 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
240 fEndOuterRadius[0] : fEndOuterRadius[1]);
241 G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
242 fEndInnerRadius[0] : fEndInnerRadius[1]);
243 G4double maxphi = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
244 std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
245
246 //
247 // Choose phi size of our segment(s) based on constants as
248 // defined in meshdefs.hh
249 //
250 // G4int numPhi = kMaxMeshSections;
251 G4double sigPhi = 2*maxphi + fDPhi;
252 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
253 G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
254
255 //
256 // We work around in phi building polygons along the way.
257 // As a reasonable compromise between accuracy and
258 // complexity (=cpu time), the following facets are chosen:
259 //
260 // 1. If fOuterRadius/maxEndOuterRad > 0.95, approximate
261 // the outer surface as a cylinder, and use one
262 // rectangular polygon (0-1) to build its mesh.
263 //
264 // Otherwise, use two trapazoidal polygons that
265 // meet at z = 0 (0-4-1)
266 //
267 // 2. If there is no inner surface, then use one
268 // polygon for each entire endcap. (0) and (1)
269 //
270 // Otherwise, use a trapazoidal polygon for each
271 // phi segment of each endcap. (0-2) and (1-3)
272 //
273 // 3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95,
274 // approximate the inner surface as a cylinder of
275 // radius fInnerRadius and use one rectangular polygon
276 // to build each phi segment of its mesh. (2-3)
277 //
278 // Otherwise, use one rectangular polygon centered
279 // at z = 0 (5-6) and two connecting trapazoidal polygons
280 // for each phi segment (2-5) and (3-6).
281 //
282
283 G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
284 G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
285
286 //
287 // Vertex assignments (v and w arrays)
288 // [0] and [1] are mandatory
289 // the rest are optional
290 //
291 // + -
292 // [0]------[4]------[1] <--- outer radius
293 // | |
294 // | |
295 // [2]---[5]---[6]---[3] <--- inner radius
296 //
297
298 G4ClippablePolygon endPoly1, endPoly2;
299
300 G4double phimax = maxphi + 0.5*fDPhi;
301 G4double phimin = - phimax;
302
303 G4ThreeVector v0, v1, v2, v3, v4, v5, v6; // -ve phi verticies for polygon
304 G4ThreeVector w0, w1, w2, w3, w4, w5, w6; // +ve phi verticies for polygon
305
306 //
307 // decide verticies of -ve phi boundary
308 //
309
310 G4double cosPhi = std::cos(phimin);
311 G4double sinPhi = std::sin(phimin);
312
313 // Outer hyperbolic surface
314
315 v0 = transform.TransformPoint(
316 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
317 + fZHalfLength));
318 v1 = transform.TransformPoint(
319 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
320 - fZHalfLength));
321 if (splitOuter)
322 {
323 v4 = transform.TransformPoint(
324 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
325 }
326
327 // Inner hyperbolic surface
328
329 G4double zInnerSplit = 0.;
330 if (splitInner)
331 {
332 v2 = transform.TransformPoint(
333 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
334 + fZHalfLength));
335 v3 = transform.TransformPoint(
336 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
337 - fZHalfLength));
338
339 // Find intersection of tangential line of inner
340 // surface at z = fZHalfLength and line r=fInnerRadius.
341 G4double dr = fZHalfLength * fTanInnerStereo2;
342 G4double dz = maxEndInnerRad;
343 zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
344
345 // Build associated vertices
346 v5 = transform.TransformPoint(
347 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
348 + zInnerSplit));
349 v6 = transform.TransformPoint(
350 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
351 - zInnerSplit));
352 }
353 else
354 {
355 v2 = transform.TransformPoint(
356 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
357 + fZHalfLength));
358 v3 = transform.TransformPoint(
359 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
360 - fZHalfLength));
361 }
362
363 //
364 // decide vertices of +ve phi boundary
365 //
366
367 cosPhi = std::cos(phimax);
368 sinPhi = std::sin(phimax);
369
370 // Outer hyperbolic surface
371
372 w0 = transform.TransformPoint(
373 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
374 + fZHalfLength));
375 w1 = transform.TransformPoint(
376 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
377 - fZHalfLength));
378 if (splitOuter)
379 {
380 G4double r = rFudge*fOuterRadius;
381
382 w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
383
384 AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
385 AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
386 }
387 else
388 {
389 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
390 }
391
392 // Inner hyperbolic surface
393
394 if (splitInner)
395 {
396 w2 = transform.TransformPoint(
397 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
398 + fZHalfLength));
399 w3 = transform.TransformPoint(
400 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
401 - fZHalfLength));
402
403 w5 = transform.TransformPoint(
404 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
405 + zInnerSplit));
406 w6 = transform.TransformPoint(
407 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
408 - zInnerSplit));
409
410 AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
411 AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
412 AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
413
414 }
415 else
416 {
417 w2 = transform.TransformPoint(
418 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
419 + fZHalfLength));
420 w3 = transform.TransformPoint(
421 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
422 - fZHalfLength));
423
424 AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
425 }
426
427 //
428 // Endplate segments
429 //
430 AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
431 AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
432
433 //
434 // Return min/max value
435 //
436 return extentList.GetExtent( min, max );
437}
438
439
440//=====================================================================
441//* AddPolyToExtent ---------------------------------------------------
442
443void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
444 const G4ThreeVector &v1,
445 const G4ThreeVector &w1,
446 const G4ThreeVector &w0,
447 const G4VoxelLimits &voxelLimit,
448 const EAxis axis,
449 G4SolidExtentList &extentList )
450{
451 // Utility function for CalculateExtent
452 //
453 G4ClippablePolygon phiPoly;
454
455 phiPoly.AddVertexInOrder( v0 );
456 phiPoly.AddVertexInOrder( v1 );
457 phiPoly.AddVertexInOrder( w1 );
458 phiPoly.AddVertexInOrder( w0 );
459
460 if (phiPoly.PartialClip( voxelLimit, axis ))
461 {
462 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
463 extentList.AddSurface( phiPoly );
464 }
465}
466
467
468//=====================================================================
469//* Inside ------------------------------------------------------------
470
471EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
472{
473
474 static const G4double halftol
475 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
476 // static G4int timerid = -1;
477 // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
478 // timer.Start();
479
480
481
482 G4ThreeVector *tmpp;
483 EInside *tmpinside;
484 if (fLastInside.p == p)
485 {
486 return fLastInside.inside;
487 }
488 else
489 {
490 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
491 tmpinside = const_cast<EInside*>(&(fLastInside.inside));
492 tmpp->set(p.x(), p.y(), p.z());
493 }
494
495 EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
496 G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
497 G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
498
499 if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
500 {
501 *tmpinside = kOutside;
502 }
503 else if (outerhypearea == kSurface)
504 {
505 *tmpinside = kSurface;
506 }
507 else
508 {
509 if (distanceToOut <= halftol)
510 {
511 *tmpinside = kSurface;
512 }
513 else
514 {
515 *tmpinside = kInside;
516 }
517 }
518
519 return fLastInside.inside;
520}
521
522//=====================================================================
523//* SurfaceNormal -----------------------------------------------------
524
525G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
526{
527 //
528 // return the normal unit vector to the Hyperbolical Surface at a point
529 // p on (or nearly on) the surface
530 //
531 // Which of the three or four surfaces are we closest to?
532 //
533
534 if (fLastNormal.p == p)
535 {
536 return fLastNormal.vec;
537 }
538 G4ThreeVector *tmpp =
539 const_cast<G4ThreeVector*>(&(fLastNormal.p));
540 G4ThreeVector *tmpnormal =
541 const_cast<G4ThreeVector*>(&(fLastNormal.vec));
542 G4VTwistSurface **tmpsurface =
543 const_cast<G4VTwistSurface**>(fLastNormal.surface);
544 tmpp->set(p.x(), p.y(), p.z());
545
546 G4double distance = kInfinity;
547
548 G4VTwistSurface *surfaces[6];
549 surfaces[0] = fLatterTwisted;
550 surfaces[1] = fFormerTwisted;
551 surfaces[2] = fInnerHype;
552 surfaces[3] = fOuterHype;
553 surfaces[4] = fLowerEndcap;
554 surfaces[5] = fUpperEndcap;
555
556 G4ThreeVector xx;
557 G4ThreeVector bestxx;
558 G4int i;
559 G4int besti = -1;
560 for (i=0; i< 6; i++)
561 {
562 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
563 if (tmpdistance < distance)
564 {
565 distance = tmpdistance;
566 bestxx = xx;
567 besti = i;
568 }
569 }
570
571 tmpsurface[0] = surfaces[besti];
572 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
573
574 return fLastNormal.vec;
575}
576
577//=====================================================================
578//* DistanceToIn (p, v) -----------------------------------------------
579
580G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
581 const G4ThreeVector& v ) const
582{
583
584 // DistanceToIn (p, v):
585 // Calculate distance to surface of shape from `outside'
586 // along with the v, allowing for tolerance.
587 // The function returns kInfinity if no intersection or
588 // just grazing within tolerance.
589
590 //
591 // checking last value
592 //
593
594 G4ThreeVector *tmpp;
595 G4ThreeVector *tmpv;
596 G4double *tmpdist;
597 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
598 {
599 return fLastDistanceToIn.value;
600 }
601 else
602 {
603 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
604 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
605 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
606 tmpp->set(p.x(), p.y(), p.z());
607 tmpv->set(v.x(), v.y(), v.z());
608 }
609
610 //
611 // Calculate DistanceToIn(p,v)
612 //
613
614 EInside currentside = Inside(p);
615
616 if (currentside == kInside)
617 {
618 }
619 else
620 {
621 if (currentside == kSurface)
622 {
623 // particle is just on a boundary.
624 // If the particle is entering to the volume, return 0.
625 //
626 G4ThreeVector normal = SurfaceNormal(p);
627 if (normal*v < 0)
628 {
629 *tmpdist = 0;
630 return fLastDistanceToInWithV.value;
631 }
632 }
633 }
634
635 // now, we can take smallest positive distance.
636
637 // Initialize
638 //
639 G4double distance = kInfinity;
640
641 // find intersections and choose nearest one.
642 //
643 G4VTwistSurface *surfaces[6];
644 surfaces[0] = fLowerEndcap;
645 surfaces[1] = fUpperEndcap;
646 surfaces[2] = fLatterTwisted;
647 surfaces[3] = fFormerTwisted;
648 surfaces[4] = fInnerHype;
649 surfaces[5] = fOuterHype;
650
651 G4ThreeVector xx;
652 G4ThreeVector bestxx;
653 G4int i;
654 G4int besti = -1;
655 for (i=0; i< 6; i++)
656 {
657 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
658 if (tmpdistance < distance)
659 {
660 distance = tmpdistance;
661 bestxx = xx;
662 besti = i;
663 }
664 }
665 *tmpdist = distance;
666
667 return fLastDistanceToInWithV.value;
668}
669
670//=====================================================================
671//* DistanceToIn (p) --------------------------------------------------
672
673G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
674{
675 // DistanceToIn(p):
676 // Calculate distance to surface of shape from `outside',
677 // allowing for tolerance
678
679 //
680 // checking last value
681 //
682
683 G4ThreeVector *tmpp;
684 G4double *tmpdist;
685 if (fLastDistanceToIn.p == p)
686 {
687 return fLastDistanceToIn.value;
688 }
689 else
690 {
691 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
692 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
693 tmpp->set(p.x(), p.y(), p.z());
694 }
695
696 //
697 // Calculate DistanceToIn(p)
698 //
699
700 EInside currentside = Inside(p);
701
702 switch (currentside)
703 {
704 case (kInside) :
705 {}
706 case (kSurface) :
707 {
708 *tmpdist = 0.;
709 return fLastDistanceToIn.value;
710 }
711 case (kOutside) :
712 {
713 // Initialize
714 G4double distance = kInfinity;
715
716 // find intersections and choose nearest one.
717 G4VTwistSurface *surfaces[6];
718 surfaces[0] = fLowerEndcap;
719 surfaces[1] = fUpperEndcap;
720 surfaces[2] = fLatterTwisted;
721 surfaces[3] = fFormerTwisted;
722 surfaces[4] = fInnerHype;
723 surfaces[5] = fOuterHype;
724
725 G4int i;
726 G4int besti = -1;
727 G4ThreeVector xx;
728 G4ThreeVector bestxx;
729 for (i=0; i< 6; i++)
730 {
731 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
732 if (tmpdistance < distance)
733 {
734 distance = tmpdistance;
735 bestxx = xx;
736 besti = i;
737 }
738 }
739 *tmpdist = distance;
740 return fLastDistanceToIn.value;
741 }
742 default :
743 {
744 G4Exception("G4TwistedTubs::DistanceToIn(p)", "InvalidCondition",
745 FatalException, "Unknown point location!");
746 }
747 } // switch end
748
749 return kInfinity;
750}
751
752//=====================================================================
753//* DistanceToOut (p, v) ----------------------------------------------
754
755G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
756 const G4ThreeVector& v,
757 const G4bool calcNorm,
758 G4bool *validNorm,
759 G4ThreeVector *norm ) const
760{
761 // DistanceToOut (p, v):
762 // Calculate distance to surface of shape from `inside'
763 // along with the v, allowing for tolerance.
764 // The function returns kInfinity if no intersection or
765 // just grazing within tolerance.
766
767 //
768 // checking last value
769 //
770
771 G4ThreeVector *tmpp;
772 G4ThreeVector *tmpv;
773 G4double *tmpdist;
774 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
775 {
776 return fLastDistanceToOutWithV.value;
777 }
778 else
779 {
780 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
781 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
782 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
783 tmpp->set(p.x(), p.y(), p.z());
784 tmpv->set(v.x(), v.y(), v.z());
785 }
786
787 //
788 // Calculate DistanceToOut(p,v)
789 //
790
791 EInside currentside = Inside(p);
792
793 if (currentside == kOutside)
794 {
795 }
796 else
797 {
798 if (currentside == kSurface)
799 {
800 // particle is just on a boundary.
801 // If the particle is exiting from the volume, return 0.
802 //
803 G4ThreeVector normal = SurfaceNormal(p);
804 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
805 if (normal*v > 0)
806 {
807 if (calcNorm)
808 {
809 *norm = (blockedsurface->GetNormal(p, true));
810 *validNorm = blockedsurface->IsValidNorm();
811 }
812 *tmpdist = 0.;
813 return fLastDistanceToOutWithV.value;
814 }
815 }
816 }
817
818 // now, we can take smallest positive distance.
819
820 // Initialize
821 //
822 G4double distance = kInfinity;
823
824 // find intersections and choose nearest one.
825 //
826 G4VTwistSurface *surfaces[6];
827 surfaces[0] = fLatterTwisted;
828 surfaces[1] = fFormerTwisted;
829 surfaces[2] = fInnerHype;
830 surfaces[3] = fOuterHype;
831 surfaces[4] = fLowerEndcap;
832 surfaces[5] = fUpperEndcap;
833
834 G4int i;
835 G4int besti = -1;
836 G4ThreeVector xx;
837 G4ThreeVector bestxx;
838 for (i=0; i< 6; i++)
839 {
840 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
841 if (tmpdistance < distance)
842 {
843 distance = tmpdistance;
844 bestxx = xx;
845 besti = i;
846 }
847 }
848
849 if (calcNorm)
850 {
851 if (besti != -1)
852 {
853 *norm = (surfaces[besti]->GetNormal(p, true));
854 *validNorm = surfaces[besti]->IsValidNorm();
855 }
856 }
857
858 *tmpdist = distance;
859
860 return fLastDistanceToOutWithV.value;
861}
862
863
864//=====================================================================
865//* DistanceToOut (p) ----------------------------------------------
866
867G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
868{
869 // DistanceToOut(p):
870 // Calculate distance to surface of shape from `inside',
871 // allowing for tolerance
872
873 //
874 // checking last value
875 //
876
877 G4ThreeVector *tmpp;
878 G4double *tmpdist;
879 if (fLastDistanceToOut.p == p)
880 {
881 return fLastDistanceToOut.value;
882 }
883 else
884 {
885 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
886 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
887 tmpp->set(p.x(), p.y(), p.z());
888 }
889
890 //
891 // Calculate DistanceToOut(p)
892 //
893
894 EInside currentside = Inside(p);
895
896 switch (currentside)
897 {
898 case (kOutside) :
899 {
900 }
901 case (kSurface) :
902 {
903 *tmpdist = 0.;
904 return fLastDistanceToOut.value;
905 }
906 case (kInside) :
907 {
908 // Initialize
909 G4double distance = kInfinity;
910
911 // find intersections and choose nearest one.
912 G4VTwistSurface *surfaces[6];
913 surfaces[0] = fLatterTwisted;
914 surfaces[1] = fFormerTwisted;
915 surfaces[2] = fInnerHype;
916 surfaces[3] = fOuterHype;
917 surfaces[4] = fLowerEndcap;
918 surfaces[5] = fUpperEndcap;
919
920 G4int i;
921 G4int besti = -1;
922 G4ThreeVector xx;
923 G4ThreeVector bestxx;
924 for (i=0; i< 6; i++)
925 {
926 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
927 if (tmpdistance < distance)
928 {
929 distance = tmpdistance;
930 bestxx = xx;
931 besti = i;
932 }
933 }
934 *tmpdist = distance;
935
936 return fLastDistanceToOut.value;
937 }
938 default :
939 {
940 G4Exception("G4TwistedTubs::DistanceToOut(p)", "InvalidCondition",
941 FatalException, "Unknown point location!");
942 }
943 } // switch end
944
945 return 0;
946}
947
948//=====================================================================
949//* StreamInfo --------------------------------------------------------
950
951std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
952{
953 //
954 // Stream object contents to an output stream
955 //
956 os << "-----------------------------------------------------------\n"
957 << " *** Dump for solid - " << GetName() << " ***\n"
958 << " ===================================================\n"
959 << " Solid type: G4TwistedTubs\n"
960 << " Parameters: \n"
961 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
962 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
963 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
964 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
965 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
966 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
967 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
968 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
969 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
970 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
971 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
972 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
973 << "-----------------------------------------------------------\n";
974
975 return os;
976}
977
978
979//=====================================================================
980//* DiscribeYourselfTo ------------------------------------------------
981
982void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
983{
984 scene.AddSolid (*this);
985}
986
987//=====================================================================
988//* GetExtent ---------------------------------------------------------
989
990G4VisExtent G4TwistedTubs::GetExtent() const
991{
992 // Define the sides of the box into which the G4Tubs instance would fit.
993
994 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
995 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
996 -maxEndOuterRad, maxEndOuterRad,
997 -fZHalfLength, fZHalfLength );
998}
999
1000//=====================================================================
1001//* CreatePolyhedron --------------------------------------------------
1002
1003G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
1004{
1005 // number of meshes
1006 //
1007 G4double dA = std::max(fDPhi,fPhiTwist);
1008 const G4int m =
1009 G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
1010 const G4int n =
1011 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1012
1013 const G4int nnodes = 4*(m-1)*(n-2) + 2*m*m ;
1014 const G4int nfaces = 4*(m-1)*(n-1) + 2*(m-1)*(m-1) ;
1015
1016 G4Polyhedron *ph=new G4Polyhedron;
1017 typedef G4double G4double3[3];
1018 typedef G4int G4int4[4];
1019 G4double3* xyz = new G4double3[nnodes]; // number of nodes
1020 G4int4* faces = new G4int4[nfaces] ; // number of faces
1021 fLowerEndcap->GetFacets(m,m,xyz,faces,0) ;
1022 fUpperEndcap->GetFacets(m,m,xyz,faces,1) ;
1023 fInnerHype->GetFacets(m,n,xyz,faces,2) ;
1024 fFormerTwisted->GetFacets(m,n,xyz,faces,3) ;
1025 fOuterHype->GetFacets(m,n,xyz,faces,4) ;
1026 fLatterTwisted->GetFacets(m,n,xyz,faces,5) ;
1027
1028 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1029
1030 delete[] xyz;
1031 delete[] faces;
1032
1033 return ph;
1034}
1035
1036//=====================================================================
1037//* CreateNUBS --------------------------------------------------------
1038
1039G4NURBS* G4TwistedTubs::CreateNURBS () const
1040{
1041 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1042 G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1043 return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength);
1044 // Tube for now!!!
1045}
1046
1047//=====================================================================
1048//* GetPolyhedron -----------------------------------------------------
1049
1050G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
1051{
1052 if ((!fpPolyhedron) ||
1053 (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1054 fpPolyhedron->GetNumberOfRotationSteps()))
1055 {
1056 delete fpPolyhedron;
1057 fpPolyhedron = CreatePolyhedron();
1058 }
1059 return fpPolyhedron;
1060}
1061
1062//=====================================================================
1063//* CreateSurfaces ----------------------------------------------------
1064
1065void G4TwistedTubs::CreateSurfaces()
1066{
1067 // create 6 surfaces of TwistedTub
1068
1069 G4ThreeVector x0(0, 0, fEndZ[0]);
1070 G4ThreeVector n (0, 0, -1);
1071
1072 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
1073 fEndInnerRadius, fEndOuterRadius,
1074 fDPhi, fEndPhi, fEndZ, -1) ;
1075
1076 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
1077 fEndInnerRadius, fEndOuterRadius,
1078 fDPhi, fEndPhi, fEndZ, 1) ;
1079
1080 G4RotationMatrix rotHalfDPhi;
1081 rotHalfDPhi.rotateZ(0.5*fDPhi);
1082
1083 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
1084 fEndInnerRadius, fEndOuterRadius,
1085 fDPhi, fEndPhi, fEndZ,
1086 fInnerRadius, fOuterRadius, fKappa,
1087 1 ) ;
1088 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
1089 fEndInnerRadius, fEndOuterRadius,
1090 fDPhi, fEndPhi, fEndZ,
1091 fInnerRadius, fOuterRadius, fKappa,
1092 -1 ) ;
1093
1094 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
1095 fEndInnerRadius, fEndOuterRadius,
1096 fDPhi, fEndPhi, fEndZ,
1097 fInnerRadius, fOuterRadius,fKappa,
1098 fTanInnerStereo, fTanOuterStereo, -1) ;
1099 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
1100 fEndInnerRadius, fEndOuterRadius,
1101 fDPhi, fEndPhi, fEndZ,
1102 fInnerRadius, fOuterRadius,fKappa,
1103 fTanInnerStereo, fTanOuterStereo, 1) ;
1104
1105
1106 // set neighbour surfaces
1107 //
1108 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1109 fOuterHype, fFormerTwisted);
1110 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1111 fOuterHype, fFormerTwisted);
1112 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1113 fOuterHype, fUpperEndcap);
1114 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1115 fOuterHype, fUpperEndcap);
1116 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1117 fFormerTwisted, fUpperEndcap);
1118 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1119 fFormerTwisted, fUpperEndcap);
1120}
1121
1122
1123//=====================================================================
1124//* GetEntityType -----------------------------------------------------
1125
1126G4GeometryType G4TwistedTubs::GetEntityType() const
1127{
1128 return G4String("G4TwistedTubs");
1129}
1130
1131//=====================================================================
1132//* GetCubicVolume ----------------------------------------------------
1133
1134G4double G4TwistedTubs::GetCubicVolume()
1135{
1136 if(fCubicVolume != 0.) {;}
1137 else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1138 -fInnerRadius*fInnerRadius); }
1139 return fCubicVolume;
1140}
1141
1142//=====================================================================
1143//* GetSurfaceArea ----------------------------------------------------
1144
1145G4double G4TwistedTubs::GetSurfaceArea()
1146{
1147 if(fSurfaceArea != 0.) {;}
1148 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1149 return fSurfaceArea;
1150}
1151
1152//=====================================================================
1153//* GetPointOnSurface -------------------------------------------------
1154
1155G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1156{
1157
1158 G4double z = CLHEP::RandFlat::shoot(fEndZ[0],fEndZ[1]);
1159 G4double phi , phimin, phimax ;
1160 G4double x , xmin, xmax ;
1161 G4double r , rmin, rmax ;
1162
1163 G4double a1 = fOuterHype->GetSurfaceArea() ;
1164 G4double a2 = fInnerHype->GetSurfaceArea() ;
1165 G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1166 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1167 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1168 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1169
1170 G4double chose = CLHEP::RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1171
1172 if(chose < a1)
1173 {
1174
1175 phimin = fOuterHype->GetBoundaryMin(z) ;
1176 phimax = fOuterHype->GetBoundaryMax(z) ;
1177 phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1178
1179 return fOuterHype->SurfacePoint(phi,z,true) ;
1180
1181 }
1182 else if ( (chose >= a1) && (chose < a1 + a2 ) )
1183 {
1184
1185 phimin = fInnerHype->GetBoundaryMin(z) ;
1186 phimax = fInnerHype->GetBoundaryMax(z) ;
1187 phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1188
1189 return fInnerHype->SurfacePoint(phi,z,true) ;
1190
1191 }
1192 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1193 {
1194
1195 xmin = fLatterTwisted->GetBoundaryMin(z) ;
1196 xmax = fLatterTwisted->GetBoundaryMax(z) ;
1197 x = CLHEP::RandFlat::shoot(xmin,xmax) ;
1198
1199 return fLatterTwisted->SurfacePoint(x,z,true) ;
1200
1201 }
1202 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1203 {
1204
1205 xmin = fFormerTwisted->GetBoundaryMin(z) ;
1206 xmax = fFormerTwisted->GetBoundaryMax(z) ;
1207 x = CLHEP::RandFlat::shoot(xmin,xmax) ;
1208
1209 return fFormerTwisted->SurfacePoint(x,z,true) ;
1210
1211 }
1212 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1213 {
1214
1215 rmin = GetEndInnerRadius(0) ;
1216 rmax = GetEndOuterRadius(0) ;
1217 r = CLHEP::RandFlat::shoot(rmin,rmax) ;
1218
1219 phimin = fLowerEndcap->GetBoundaryMin(r) ;
1220 phimax = fLowerEndcap->GetBoundaryMax(r) ;
1221 phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1222
1223 return fLowerEndcap->SurfacePoint(phi,r,true) ;
1224
1225 }
1226 else
1227 {
1228 rmin = GetEndInnerRadius(1) ;
1229 rmax = GetEndOuterRadius(1) ;
1230 r = CLHEP::RandFlat::shoot(rmin,rmax) ;
1231
1232 phimin = fUpperEndcap->GetBoundaryMin(r) ;
1233 phimax = fUpperEndcap->GetBoundaryMax(r) ;
1234 phi = CLHEP::RandFlat::shoot(phimin,phimax) ;
1235
1236 return fUpperEndcap->SurfacePoint(phi,r,true) ;
1237 }
1238}
Note: See TracBrowser for help on using the repository browser.