Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TwistedTubs.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4TwistedTubs implementation
27//
28// 01-Aug-2002 - Kotoyo Hoshina ([email protected]), created.
29// 13-Nov-2003 - O.Link ([email protected]), Integration in Geant4
30// from original version in Jupiter-2.5.02 application.
31// --------------------------------------------------------------------
32
33#include "G4TwistedTubs.hh"
34
36#include "G4SystemOfUnits.hh"
38#include "G4VoxelLimits.hh"
39#include "G4AffineTransform.hh"
40#include "G4BoundingEnvelope.hh"
41#include "G4ClippablePolygon.hh"
43#include "meshdefs.hh"
44
45#include "G4VGraphicsScene.hh"
46#include "G4Polyhedron.hh"
47#include "G4VisExtent.hh"
48
49#include "Randomize.hh"
50
51#include "G4AutoLock.hh"
52
53namespace
54{
55 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
56}
57
58//=====================================================================
59//* constructors ------------------------------------------------------
60
62 G4double twistedangle,
63 G4double endinnerrad,
64 G4double endouterrad,
65 G4double halfzlen,
66 G4double dphi)
67 : G4VSolid(pname), fDPhi(dphi),
68 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
69 fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
70{
71 if (endinnerrad < DBL_MIN)
72 {
73 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
74 FatalErrorInArgument, "Invalid end-inner-radius!");
75 }
76
77 G4double sinhalftwist = std::sin(0.5 * twistedangle);
78
79 G4double endinnerradX = endinnerrad * sinhalftwist;
80 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
81 - endinnerradX * endinnerradX );
82
83 G4double endouterradX = endouterrad * sinhalftwist;
84 G4double outerrad = std::sqrt( endouterrad * endouterrad
85 - endouterradX * endouterradX );
86
87 // temporary treatment!!
88 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
89 CreateSurfaces();
90}
91
93 G4double twistedangle,
94 G4double endinnerrad,
95 G4double endouterrad,
96 G4double halfzlen,
97 G4int nseg,
98 G4double totphi)
99 : G4VSolid(pname),
100 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
101 fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
102{
103
104 if (!nseg)
105 {
106 std::ostringstream message;
107 message << "Invalid number of segments." << G4endl
108 << " nseg = " << nseg;
109 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
110 FatalErrorInArgument, message);
111 }
112 if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
113 {
114 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
115 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
116 }
117
118 G4double sinhalftwist = std::sin(0.5 * twistedangle);
119
120 G4double endinnerradX = endinnerrad * sinhalftwist;
121 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
122 - endinnerradX * endinnerradX );
123
124 G4double endouterradX = endouterrad * sinhalftwist;
125 G4double outerrad = std::sqrt( endouterrad * endouterrad
126 - endouterradX * endouterradX );
127
128 // temporary treatment!!
129 fDPhi = totphi / nseg;
130 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
131 CreateSurfaces();
132}
133
135 G4double twistedangle,
136 G4double innerrad,
137 G4double outerrad,
138 G4double negativeEndz,
139 G4double positiveEndz,
140 G4double dphi)
141 : G4VSolid(pname), fDPhi(dphi),
142 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
143 fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
144{
145 if (innerrad < DBL_MIN)
146 {
147 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
148 FatalErrorInArgument, "Invalid end-inner-radius!");
149 }
150
151 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
152 CreateSurfaces();
153}
154
156 G4double twistedangle,
157 G4double innerrad,
158 G4double outerrad,
159 G4double negativeEndz,
160 G4double positiveEndz,
161 G4int nseg,
162 G4double totphi)
163 : G4VSolid(pname),
164 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
165 fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
166{
167 if (!nseg)
168 {
169 std::ostringstream message;
170 message << "Invalid number of segments." << G4endl
171 << " nseg = " << nseg;
172 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
173 FatalErrorInArgument, message);
174 }
175 if (totphi == DBL_MIN || innerrad < DBL_MIN)
176 {
177 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
178 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
179 }
180
181 fDPhi = totphi / nseg;
182 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
183 CreateSurfaces();
184}
185
186//=====================================================================
187//* Fake default constructor ------------------------------------------
188
190 : G4VSolid(a), fPhiTwist(0.), fInnerRadius(0.), fOuterRadius(0.), fDPhi(0.),
191 fZHalfLength(0.), fInnerStereo(0.), fOuterStereo(0.), fTanInnerStereo(0.),
192 fTanOuterStereo(0.), fKappa(0.), fInnerRadius2(0.), fOuterRadius2(0.),
193 fTanInnerStereo2(0.), fTanOuterStereo2(0.), fLowerEndcap(0), fUpperEndcap(0),
194 fLatterTwisted(0), fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
195{
196}
197
198//=====================================================================
199//* destructor --------------------------------------------------------
200
202{
203 if (fLowerEndcap) { delete fLowerEndcap; }
204 if (fUpperEndcap) { delete fUpperEndcap; }
205 if (fLatterTwisted) { delete fLatterTwisted; }
206 if (fFormerTwisted) { delete fFormerTwisted; }
207 if (fInnerHype) { delete fInnerHype; }
208 if (fOuterHype) { delete fOuterHype; }
209 if (fpPolyhedron) { delete fpPolyhedron; fpPolyhedron = nullptr; }
210}
211
212//=====================================================================
213//* Copy constructor --------------------------------------------------
214
216 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
217 fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
218 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
219 fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
220 fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
221 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
222 fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
223 fTanOuterStereo2(rhs.fTanOuterStereo2),
224 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), fFormerTwisted(0),
225 fInnerHype(0), fOuterHype(0),
226 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
227 fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
228 fLastDistanceToIn(rhs.fLastDistanceToIn),
229 fLastDistanceToOut(rhs.fLastDistanceToOut),
230 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
231 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
232{
233 for (auto i=0; i<2; ++i)
234 {
235 fEndZ[i] = rhs.fEndZ[i];
236 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
237 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
238 fEndPhi[i] = rhs.fEndPhi[i];
239 fEndZ2[i] = rhs.fEndZ2[i];
240 }
241 CreateSurfaces();
242}
243
244
245//=====================================================================
246//* Assignment operator -----------------------------------------------
247
249{
250 // Check assignment to self
251 //
252 if (this == &rhs) { return *this; }
253
254 // Copy base class data
255 //
257
258 // Copy data
259 //
260 fPhiTwist= rhs.fPhiTwist;
261 fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
262 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
263 fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
264 fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
265 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
266 fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
267 fTanOuterStereo2= rhs.fTanOuterStereo2;
268 fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= 0;
269 fInnerHype= fOuterHype= 0;
270 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
271 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
272 fLastDistanceToIn= rhs.fLastDistanceToIn;
273 fLastDistanceToOut= rhs.fLastDistanceToOut;
274 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
275 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
276
277 for (auto i=0; i<2; ++i)
278 {
279 fEndZ[i] = rhs.fEndZ[i];
280 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
281 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
282 fEndPhi[i] = rhs.fEndPhi[i];
283 fEndZ2[i] = rhs.fEndZ2[i];
284 }
285
286 CreateSurfaces();
287 fRebuildPolyhedron = false;
288 delete fpPolyhedron; fpPolyhedron = nullptr;
289
290 return *this;
291}
292
293//=====================================================================
294//* ComputeDimensions -------------------------------------------------
295
297 const G4int /* n */ ,
298 const G4VPhysicalVolume* /* pRep */ )
299{
300 G4Exception("G4TwistedTubs::ComputeDimensions()",
301 "GeomSolids0001", FatalException,
302 "G4TwistedTubs does not support Parameterisation.");
303}
304
305//=====================================================================
306//* BoundingLimits ----------------------------------------------------
307
309 G4ThreeVector& pMax) const
310{
311 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
312 fEndOuterRadius[0] : fEndOuterRadius[1]);
313 pMin.set(-maxEndOuterRad,-maxEndOuterRad,-fZHalfLength);
314 pMax.set( maxEndOuterRad, maxEndOuterRad, fZHalfLength);
315
316 // Check correctness of the bounding box
317 //
318 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
319 {
320 std::ostringstream message;
321 message << "Bad bounding box (min >= max) for solid: "
322 << GetName() << " !"
323 << "\npMin = " << pMin
324 << "\npMax = " << pMax;
325 G4Exception("G4TwistedTubs::BoundingLimits()", "GeomMgt0001",
326 JustWarning, message);
327 DumpInfo();
328 }
329}
330
331//=====================================================================
332//* CalculateExtent ---------------------------------------------------
333
334G4bool
336 const G4VoxelLimits& pVoxelLimit,
337 const G4AffineTransform& pTransform,
338 G4double& pMin, G4double& pMax) const
339{
340 G4ThreeVector bmin, bmax;
341
342 // Get bounding box
343 BoundingLimits(bmin,bmax);
344
345 // Find extent
346 G4BoundingEnvelope bbox(bmin,bmax);
347 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
348}
349
350
351//=====================================================================
352//* Inside ------------------------------------------------------------
353
355{
356
357 const G4double halftol
359 // static G4int timerid = -1;
360 // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
361 // timer.Start();
362
363 G4ThreeVector *tmpp;
364 EInside *tmpinside;
365 if (fLastInside.p == p)
366 {
367 return fLastInside.inside;
368 }
369 else
370 {
371 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
372 tmpinside = const_cast<EInside*>(&(fLastInside.inside));
373 tmpp->set(p.x(), p.y(), p.z());
374 }
375
376 EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
377 G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
378 G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
379
380 if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
381 {
382 *tmpinside = kOutside;
383 }
384 else if (outerhypearea == kSurface)
385 {
386 *tmpinside = kSurface;
387 }
388 else
389 {
390 if (distanceToOut <= halftol)
391 {
392 *tmpinside = kSurface;
393 }
394 else
395 {
396 *tmpinside = kInside;
397 }
398 }
399
400 return fLastInside.inside;
401}
402
403//=====================================================================
404//* SurfaceNormal -----------------------------------------------------
405
407{
408 //
409 // return the normal unit vector to the Hyperbolical Surface at a point
410 // p on (or nearly on) the surface
411 //
412 // Which of the three or four surfaces are we closest to?
413 //
414
415 if (fLastNormal.p == p)
416 {
417 return fLastNormal.vec;
418 }
419 G4ThreeVector *tmpp =
420 const_cast<G4ThreeVector*>(&(fLastNormal.p));
421 G4ThreeVector *tmpnormal =
422 const_cast<G4ThreeVector*>(&(fLastNormal.vec));
423 G4VTwistSurface **tmpsurface =
424 const_cast<G4VTwistSurface**>(fLastNormal.surface);
425 tmpp->set(p.x(), p.y(), p.z());
426
427 G4double distance = kInfinity;
428
429 G4VTwistSurface *surfaces[6];
430 surfaces[0] = fLatterTwisted;
431 surfaces[1] = fFormerTwisted;
432 surfaces[2] = fInnerHype;
433 surfaces[3] = fOuterHype;
434 surfaces[4] = fLowerEndcap;
435 surfaces[5] = fUpperEndcap;
436
437 G4ThreeVector xx;
438 G4ThreeVector bestxx;
439 G4int besti = -1;
440 for (auto i=0; i<6; ++i)
441 {
442 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
443 if (tmpdistance < distance)
444 {
445 distance = tmpdistance;
446 bestxx = xx;
447 besti = i;
448 }
449 }
450
451 tmpsurface[0] = surfaces[besti];
452 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
453
454 return fLastNormal.vec;
455}
456
457//=====================================================================
458//* DistanceToIn (p, v) -----------------------------------------------
459
461 const G4ThreeVector& v ) const
462{
463
464 // DistanceToIn (p, v):
465 // Calculate distance to surface of shape from `outside'
466 // along with the v, allowing for tolerance.
467 // The function returns kInfinity if no intersection or
468 // just grazing within tolerance.
469
470 //
471 // checking last value
472 //
473
474 G4ThreeVector* tmpp;
475 G4ThreeVector* tmpv;
476 G4double* tmpdist;
477 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
478 {
479 return fLastDistanceToIn.value;
480 }
481 else
482 {
483 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
484 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
485 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
486 tmpp->set(p.x(), p.y(), p.z());
487 tmpv->set(v.x(), v.y(), v.z());
488 }
489
490 //
491 // Calculate DistanceToIn(p,v)
492 //
493
494 EInside currentside = Inside(p);
495
496 if (currentside == kInside)
497 {
498 }
499 else
500 {
501 if (currentside == kSurface)
502 {
503 // particle is just on a boundary.
504 // If the particle is entering to the volume, return 0.
505 //
506 G4ThreeVector normal = SurfaceNormal(p);
507 if (normal*v < 0)
508 {
509 *tmpdist = 0.;
510 return fLastDistanceToInWithV.value;
511 }
512 }
513 }
514
515 // now, we can take smallest positive distance.
516
517 // Initialize
518 //
519 G4double distance = kInfinity;
520
521 // find intersections and choose nearest one.
522 //
523 G4VTwistSurface* surfaces[6];
524 surfaces[0] = fLowerEndcap;
525 surfaces[1] = fUpperEndcap;
526 surfaces[2] = fLatterTwisted;
527 surfaces[3] = fFormerTwisted;
528 surfaces[4] = fInnerHype;
529 surfaces[5] = fOuterHype;
530
531 G4ThreeVector xx;
532 G4ThreeVector bestxx;
533 for (auto i=0; i<6; ++i)
534 {
535 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
536 if (tmpdistance < distance)
537 {
538 distance = tmpdistance;
539 bestxx = xx;
540 }
541 }
542 *tmpdist = distance;
543
544 return fLastDistanceToInWithV.value;
545}
546
547//=====================================================================
548//* DistanceToIn (p) --------------------------------------------------
549
551{
552 // DistanceToIn(p):
553 // Calculate distance to surface of shape from `outside',
554 // allowing for tolerance
555
556 //
557 // checking last value
558 //
559
560 G4ThreeVector* tmpp;
561 G4double* tmpdist;
562 if (fLastDistanceToIn.p == p)
563 {
564 return fLastDistanceToIn.value;
565 }
566 else
567 {
568 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
569 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
570 tmpp->set(p.x(), p.y(), p.z());
571 }
572
573 //
574 // Calculate DistanceToIn(p)
575 //
576
577 EInside currentside = Inside(p);
578
579 switch (currentside)
580 {
581 case (kInside) :
582 {}
583 case (kSurface) :
584 {
585 *tmpdist = 0.;
586 return fLastDistanceToIn.value;
587 }
588 case (kOutside) :
589 {
590 // Initialize
591 G4double distance = kInfinity;
592
593 // find intersections and choose nearest one.
594 G4VTwistSurface *surfaces[6];
595 surfaces[0] = fLowerEndcap;
596 surfaces[1] = fUpperEndcap;
597 surfaces[2] = fLatterTwisted;
598 surfaces[3] = fFormerTwisted;
599 surfaces[4] = fInnerHype;
600 surfaces[5] = fOuterHype;
601
602 G4ThreeVector xx;
603 G4ThreeVector bestxx;
604 for (auto i=0; i<6; ++i)
605 {
606 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
607 if (tmpdistance < distance)
608 {
609 distance = tmpdistance;
610 bestxx = xx;
611 }
612 }
613 *tmpdist = distance;
614 return fLastDistanceToIn.value;
615 }
616 default :
617 {
618 G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
619 FatalException, "Unknown point location!");
620 }
621 } // switch end
622
623 return kInfinity;
624}
625
626//=====================================================================
627//* DistanceToOut (p, v) ----------------------------------------------
628
630 const G4ThreeVector& v,
631 const G4bool calcNorm,
632 G4bool *validNorm,
633 G4ThreeVector *norm ) const
634{
635 // DistanceToOut (p, v):
636 // Calculate distance to surface of shape from `inside'
637 // along with the v, allowing for tolerance.
638 // The function returns kInfinity if no intersection or
639 // just grazing within tolerance.
640
641 //
642 // checking last value
643 //
644
645 G4ThreeVector* tmpp;
646 G4ThreeVector* tmpv;
647 G4double* tmpdist;
648 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
649 {
650 return fLastDistanceToOutWithV.value;
651 }
652 else
653 {
654 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
655 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
656 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
657 tmpp->set(p.x(), p.y(), p.z());
658 tmpv->set(v.x(), v.y(), v.z());
659 }
660
661 //
662 // Calculate DistanceToOut(p,v)
663 //
664
665 EInside currentside = Inside(p);
666
667 if (currentside == kOutside)
668 {
669 }
670 else
671 {
672 if (currentside == kSurface)
673 {
674 // particle is just on a boundary.
675 // If the particle is exiting from the volume, return 0.
676 //
677 G4ThreeVector normal = SurfaceNormal(p);
678 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
679 if (normal*v > 0)
680 {
681 if (calcNorm)
682 {
683 *norm = (blockedsurface->GetNormal(p, true));
684 *validNorm = blockedsurface->IsValidNorm();
685 }
686 *tmpdist = 0.;
687 return fLastDistanceToOutWithV.value;
688 }
689 }
690 }
691
692 // now, we can take smallest positive distance.
693
694 // Initialize
695 //
696 G4double distance = kInfinity;
697
698 // find intersections and choose nearest one.
699 //
700 G4VTwistSurface* surfaces[6];
701 surfaces[0] = fLatterTwisted;
702 surfaces[1] = fFormerTwisted;
703 surfaces[2] = fInnerHype;
704 surfaces[3] = fOuterHype;
705 surfaces[4] = fLowerEndcap;
706 surfaces[5] = fUpperEndcap;
707
708 G4int besti = -1;
709 G4ThreeVector xx;
710 G4ThreeVector bestxx;
711 for (auto i=0; i<6; ++i)
712 {
713 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
714 if (tmpdistance < distance)
715 {
716 distance = tmpdistance;
717 bestxx = xx;
718 besti = i;
719 }
720 }
721
722 if (calcNorm)
723 {
724 if (besti != -1)
725 {
726 *norm = (surfaces[besti]->GetNormal(p, true));
727 *validNorm = surfaces[besti]->IsValidNorm();
728 }
729 }
730
731 *tmpdist = distance;
732
733 return fLastDistanceToOutWithV.value;
734}
735
736
737//=====================================================================
738//* DistanceToOut (p) ----------------------------------------------
739
741{
742 // DistanceToOut(p):
743 // Calculate distance to surface of shape from `inside',
744 // allowing for tolerance
745
746 //
747 // checking last value
748 //
749
750 G4ThreeVector* tmpp;
751 G4double* tmpdist;
752 if (fLastDistanceToOut.p == p)
753 {
754 return fLastDistanceToOut.value;
755 }
756 else
757 {
758 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
759 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
760 tmpp->set(p.x(), p.y(), p.z());
761 }
762
763 //
764 // Calculate DistanceToOut(p)
765 //
766
767 EInside currentside = Inside(p);
768
769 switch (currentside)
770 {
771 case (kOutside) :
772 {
773 }
774 case (kSurface) :
775 {
776 *tmpdist = 0.;
777 return fLastDistanceToOut.value;
778 }
779 case (kInside) :
780 {
781 // Initialize
782 G4double distance = kInfinity;
783
784 // find intersections and choose nearest one.
785 G4VTwistSurface* surfaces[6];
786 surfaces[0] = fLatterTwisted;
787 surfaces[1] = fFormerTwisted;
788 surfaces[2] = fInnerHype;
789 surfaces[3] = fOuterHype;
790 surfaces[4] = fLowerEndcap;
791 surfaces[5] = fUpperEndcap;
792
793 G4ThreeVector xx;
794 G4ThreeVector bestxx;
795 for (auto i=0; i<6; ++i)
796 {
797 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
798 if (tmpdistance < distance)
799 {
800 distance = tmpdistance;
801 bestxx = xx;
802 }
803 }
804 *tmpdist = distance;
805
806 return fLastDistanceToOut.value;
807 }
808 default :
809 {
810 G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
811 FatalException, "Unknown point location!");
812 }
813 } // switch end
814
815 return 0.;
816}
817
818//=====================================================================
819//* StreamInfo --------------------------------------------------------
820
821std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
822{
823 //
824 // Stream object contents to an output stream
825 //
826 G4int oldprc = os.precision(16);
827 os << "-----------------------------------------------------------\n"
828 << " *** Dump for solid - " << GetName() << " ***\n"
829 << " ===================================================\n"
830 << " Solid type: G4TwistedTubs\n"
831 << " Parameters: \n"
832 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
833 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
834 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
835 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
836 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
837 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
838 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
839 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
840 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
841 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
842 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
843 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
844 << "-----------------------------------------------------------\n";
845 os.precision(oldprc);
846
847 return os;
848}
849
850
851//=====================================================================
852//* DiscribeYourselfTo ------------------------------------------------
853
855{
856 scene.AddSolid (*this);
857}
858
859//=====================================================================
860//* GetExtent ---------------------------------------------------------
861
863{
864 // Define the sides of the box into which the G4Tubs instance would fit.
865
866 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
867 fEndOuterRadius[0] : fEndOuterRadius[1]);
868 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
869 -maxEndOuterRad, maxEndOuterRad,
870 -fZHalfLength, fZHalfLength );
871}
872
873//=====================================================================
874//* CreatePolyhedron --------------------------------------------------
875
877{
878 // number of meshes
879 //
880 G4double absPhiTwist = std::abs(fPhiTwist);
881 G4double dA = std::max(fDPhi,absPhiTwist);
882 const G4int k =
884 const G4int n =
885 G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
886
887 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
888 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
889
890 G4Polyhedron* ph = new G4Polyhedron;
891 typedef G4double G4double3[3];
892 typedef G4int G4int4[4];
893 G4double3* xyz = new G4double3[nnodes]; // number of nodes
894 G4int4* faces = new G4int4[nfaces] ; // number of faces
895 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
896 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
897 fInnerHype->GetFacets(k,n,xyz,faces,2) ;
898 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
899 fOuterHype->GetFacets(k,n,xyz,faces,4) ;
900 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
901
902 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
903
904 delete[] xyz;
905 delete[] faces;
906
907 return ph;
908}
909
910//=====================================================================
911//* GetPolyhedron -----------------------------------------------------
912
914{
915 if (fpPolyhedron == nullptr ||
916 fRebuildPolyhedron ||
918 fpPolyhedron->GetNumberOfRotationSteps())
919 {
920 G4AutoLock l(&polyhedronMutex);
921 delete fpPolyhedron;
922 fpPolyhedron = CreatePolyhedron();
923 fRebuildPolyhedron = false;
924 l.unlock();
925 }
926 return fpPolyhedron;
927}
928
929//=====================================================================
930//* CreateSurfaces ----------------------------------------------------
931
932void G4TwistedTubs::CreateSurfaces()
933{
934 // create 6 surfaces of TwistedTub
935
936 G4ThreeVector x0(0, 0, fEndZ[0]);
937 G4ThreeVector n (0, 0, -1);
938
939 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
940 fEndInnerRadius, fEndOuterRadius,
941 fDPhi, fEndPhi, fEndZ, -1) ;
942
943 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
944 fEndInnerRadius, fEndOuterRadius,
945 fDPhi, fEndPhi, fEndZ, 1) ;
946
947 G4RotationMatrix rotHalfDPhi;
948 rotHalfDPhi.rotateZ(0.5*fDPhi);
949
950 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
951 fEndInnerRadius, fEndOuterRadius,
952 fDPhi, fEndPhi, fEndZ,
953 fInnerRadius, fOuterRadius, fKappa,
954 1 ) ;
955 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
956 fEndInnerRadius, fEndOuterRadius,
957 fDPhi, fEndPhi, fEndZ,
958 fInnerRadius, fOuterRadius, fKappa,
959 -1 ) ;
960
961 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
962 fEndInnerRadius, fEndOuterRadius,
963 fDPhi, fEndPhi, fEndZ,
964 fInnerRadius, fOuterRadius,fKappa,
965 fTanInnerStereo, fTanOuterStereo, -1) ;
966 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
967 fEndInnerRadius, fEndOuterRadius,
968 fDPhi, fEndPhi, fEndZ,
969 fInnerRadius, fOuterRadius,fKappa,
970 fTanInnerStereo, fTanOuterStereo, 1) ;
971
972
973 // set neighbour surfaces
974 //
975 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
976 fOuterHype, fFormerTwisted);
977 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
978 fOuterHype, fFormerTwisted);
979 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
980 fOuterHype, fUpperEndcap);
981 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
982 fOuterHype, fUpperEndcap);
983 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
984 fFormerTwisted, fUpperEndcap);
985 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
986 fFormerTwisted, fUpperEndcap);
987}
988
989
990//=====================================================================
991//* GetEntityType -----------------------------------------------------
992
994{
995 return G4String("G4TwistedTubs");
996}
997
998//=====================================================================
999//* Clone -------------------------------------------------------------
1000
1002{
1003 return new G4TwistedTubs(*this);
1004}
1005
1006//=====================================================================
1007//* GetCubicVolume ----------------------------------------------------
1008
1010{
1011 if(fCubicVolume != 0.) {;}
1012 else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1013 -fInnerRadius*fInnerRadius); }
1014 return fCubicVolume;
1015}
1016
1017//=====================================================================
1018//* GetSurfaceArea ----------------------------------------------------
1019
1021{
1022 if(fSurfaceArea != 0.) {;}
1023 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1024 return fSurfaceArea;
1025}
1026
1027//=====================================================================
1028//* GetPointOnSurface -------------------------------------------------
1029
1031{
1032
1033 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1034 G4double phi , phimin, phimax ;
1035 G4double x , xmin, xmax ;
1036 G4double r , rmin, rmax ;
1037
1038 G4double a1 = fOuterHype->GetSurfaceArea() ;
1039 G4double a2 = fInnerHype->GetSurfaceArea() ;
1040 G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1041 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1042 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1043 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1044
1045 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1046
1047 if(chose < a1)
1048 {
1049
1050 phimin = fOuterHype->GetBoundaryMin(z) ;
1051 phimax = fOuterHype->GetBoundaryMax(z) ;
1052 phi = G4RandFlat::shoot(phimin,phimax) ;
1053
1054 return fOuterHype->SurfacePoint(phi,z,true) ;
1055
1056 }
1057 else if ( (chose >= a1) && (chose < a1 + a2 ) )
1058 {
1059
1060 phimin = fInnerHype->GetBoundaryMin(z) ;
1061 phimax = fInnerHype->GetBoundaryMax(z) ;
1062 phi = G4RandFlat::shoot(phimin,phimax) ;
1063
1064 return fInnerHype->SurfacePoint(phi,z,true) ;
1065
1066 }
1067 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1068 {
1069
1070 xmin = fLatterTwisted->GetBoundaryMin(z) ;
1071 xmax = fLatterTwisted->GetBoundaryMax(z) ;
1072 x = G4RandFlat::shoot(xmin,xmax) ;
1073
1074 return fLatterTwisted->SurfacePoint(x,z,true) ;
1075
1076 }
1077 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1078 {
1079
1080 xmin = fFormerTwisted->GetBoundaryMin(z) ;
1081 xmax = fFormerTwisted->GetBoundaryMax(z) ;
1082 x = G4RandFlat::shoot(xmin,xmax) ;
1083
1084 return fFormerTwisted->SurfacePoint(x,z,true) ;
1085 }
1086 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1087 {
1088 rmin = GetEndInnerRadius(0) ;
1089 rmax = GetEndOuterRadius(0) ;
1090 r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1091
1092 phimin = fLowerEndcap->GetBoundaryMin(r) ;
1093 phimax = fLowerEndcap->GetBoundaryMax(r) ;
1094 phi = G4RandFlat::shoot(phimin,phimax) ;
1095
1096 return fLowerEndcap->SurfacePoint(phi,r,true) ;
1097 }
1098 else
1099 {
1100 rmin = GetEndInnerRadius(1) ;
1101 rmax = GetEndOuterRadius(1) ;
1102 r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1103
1104 phimin = fUpperEndcap->GetBoundaryMin(r) ;
1105 phimax = fUpperEndcap->GetBoundaryMax(r) ;
1106 phi = G4RandFlat::shoot(phimin,phimax) ;
1107
1108 return fUpperEndcap->SurfacePoint(phi,r,true) ;
1109 }
1110}
@ JustWarning
@ FatalException
@ FatalErrorInArgument
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:85
std::mutex G4Mutex
Definition: G4Threading.hh:81
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
double z() const
double x() const
double y() const
double getRho() const
void set(double x, double y, double z)
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:87
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double GetRadialTolerance() const
static G4GeometryTolerance * GetInstance()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4ThreeVector GetPointOnSurface() const
G4double GetSurfaceArea()
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4Polyhedron * CreatePolyhedron() const
void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
G4VSolid * Clone() const
G4TwistedTubs & operator=(const G4TwistedTubs &rhs)
G4double GetEndInnerRadius() const
G4double GetEndOuterRadius() const
G4double GetCubicVolume()
G4Polyhedron * GetPolyhedron() const
G4GeometryType GetEntityType() const
virtual ~G4TwistedTubs()
G4TwistedTubs(const G4String &pname, G4double twistedangle, G4double endinnerrad, G4double endouterrad, G4double halfzlen, G4double dphi)
EInside Inside(const G4ThreeVector &p) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=nullptr, G4ThreeVector *n=nullptr) const
G4VisExtent GetExtent() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
std::ostream & StreamInfo(std::ostream &os) const
virtual void AddSolid(const G4Box &)=0
G4String GetName() const
void DumpInfo() const
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:98
virtual G4double GetSurfaceArea()
Definition: G4VSolid.cc:238
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
void SetNeighbours(G4VTwistSurface *ax0min, G4VTwistSurface *ax1min, G4VTwistSurface *ax0max, G4VTwistSurface *ax1max)
virtual G4double GetBoundaryMin(G4double)=0
G4bool IsValidNorm() const
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
virtual G4double GetBoundaryMax(G4double)=0
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
virtual G4double GetSurfaceArea()=0
static G4int GetNumberOfRotationSteps()
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
EAxis
Definition: geomdefs.hh:54
EInside
Definition: geomdefs.hh:67
@ kInside
Definition: geomdefs.hh:70
@ kOutside
Definition: geomdefs.hh:68
@ kSurface
Definition: geomdefs.hh:69
T sqr(const T &x)
Definition: templates.hh:128
#define DBL_MIN
Definition: templates.hh:54