| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4TwistedTubs implementation
27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch),
30 // from original version in Jupi
31 // -------------------------------------------
32
33 #include "G4TwistedTubs.hh"
34
35 #include "G4GeomTools.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4GeometryTolerance.hh"
39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42 #include "G4ClippablePolygon.hh"
43 #include "G4VPVParameterisation.hh"
44 #include "meshdefs.hh"
45
46 #include "G4VGraphicsScene.hh"
47 #include "G4Polyhedron.hh"
48 #include "G4VisExtent.hh"
49
50 #include "Randomize.hh"
51
52 #include "G4AutoLock.hh"
53
54 namespace
55 {
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
57 }
58
59
60 //============================================
61 //* constructors -----------------------------
62
63 G4TwistedTubs::G4TwistedTubs(const G4String& p
64 G4double t
65 G4double e
66 G4double e
67 G4double h
68 G4double d
69 : G4VSolid(pname), fDPhi(dphi),
70 fLowerEndcap(nullptr), fUpperEndcap(nullp
71 fFormerTwisted(nullptr), fInnerHype(nullp
72 {
73 if (endinnerrad < DBL_MIN)
74 {
75 G4Exception("G4TwistedTubs::G4TwistedTub
76 FatalErrorInArgument, "Inval
77 }
78
79 G4double sinhalftwist = std::sin(0.5 * twis
80
81 G4double endinnerradX = endinnerrad * sinha
82 G4double innerrad = std::sqrt( endinner
83 - endinnerrad
84
85 G4double endouterradX = endouterrad * sinha
86 G4double outerrad = std::sqrt( endouter
87 - endouterrad
88
89 // temporary treatment!!
90 SetFields(twistedangle, innerrad, outerrad,
91 CreateSurfaces();
92 }
93
94 G4TwistedTubs::G4TwistedTubs(const G4String& p
95 G4double t
96 G4double e
97 G4double e
98 G4double h
99 G4int n
100 G4double t
101 : G4VSolid(pname),
102 fLowerEndcap(nullptr), fUpperEndcap(nullp
103 fFormerTwisted(nullptr), fInnerHype(nullp
104 {
105
106 if (nseg == 0)
107 {
108 std::ostringstream message;
109 message << "Invalid number of segments."
110 << " nseg = " << nseg;
111 G4Exception("G4TwistedTubs::G4TwistedTub
112 FatalErrorInArgument, messag
113 }
114 if (totphi == DBL_MIN || endinnerrad < DBL_
115 {
116 G4Exception("G4TwistedTubs::G4TwistedTub
117 FatalErrorInArgument, "Invalid
118 }
119
120 G4double sinhalftwist = std::sin(0.5 * twis
121
122 G4double endinnerradX = endinnerrad * sinha
123 G4double innerrad = std::sqrt( endinner
124 - endinnerrad
125
126 G4double endouterradX = endouterrad * sinha
127 G4double outerrad = std::sqrt( endouter
128 - endouterrad
129
130 // temporary treatment!!
131 fDPhi = totphi / nseg;
132 SetFields(twistedangle, innerrad, outerrad,
133 CreateSurfaces();
134 }
135
136 G4TwistedTubs::G4TwistedTubs(const G4String& p
137 G4double t
138 G4double i
139 G4double o
140 G4double n
141 G4double p
142 G4double d
143 : G4VSolid(pname), fDPhi(dphi),
144 fLowerEndcap(nullptr), fUpperEndcap(nullp
145 fFormerTwisted(nullptr), fInnerHype(nullp
146 {
147 if (innerrad < DBL_MIN)
148 {
149 G4Exception("G4TwistedTubs::G4TwistedTub
150 FatalErrorInArgument, "Inval
151 }
152
153 SetFields(twistedangle, innerrad, outerrad,
154 CreateSurfaces();
155 }
156
157 G4TwistedTubs::G4TwistedTubs(const G4String& p
158 G4double t
159 G4double i
160 G4double o
161 G4double n
162 G4double p
163 G4int n
164 G4double t
165 : G4VSolid(pname),
166 fLowerEndcap(nullptr), fUpperEndcap(nullp
167 fFormerTwisted(nullptr), fInnerHype(nullp
168 {
169 if (nseg == 0)
170 {
171 std::ostringstream message;
172 message << "Invalid number of segments."
173 << " nseg = " << nseg;
174 G4Exception("G4TwistedTubs::G4TwistedTub
175 FatalErrorInArgument, messag
176 }
177 if (totphi == DBL_MIN || innerrad < DBL_MIN
178 {
179 G4Exception("G4TwistedTubs::G4TwistedTub
180 FatalErrorInArgument, "Invalid
181 }
182
183 fDPhi = totphi / nseg;
184 SetFields(twistedangle, innerrad, outerrad,
185 CreateSurfaces();
186 }
187
188 //============================================
189 //* Fake default constructor -----------------
190
191 G4TwistedTubs::G4TwistedTubs( __void__& a )
192 : G4VSolid(a),
193 fLowerEndcap(nullptr), fUpperEndcap(nullpt
194 fLatterTwisted(nullptr), fFormerTwisted(nu
195 fInnerHype(nullptr), fOuterHype(nullptr)
196 {
197 }
198
199 //============================================
200 //* destructor -------------------------------
201
202 G4TwistedTubs::~G4TwistedTubs()
203 {
204 delete fLowerEndcap;
205 delete fUpperEndcap;
206 delete fLatterTwisted;
207 delete fFormerTwisted;
208 delete fInnerHype;
209 delete fOuterHype;
210 delete fpPolyhedron; fpPolyhedron = nullptr
211 }
212
213 //============================================
214 //* Copy constructor -------------------------
215
216 G4TwistedTubs::G4TwistedTubs(const G4TwistedTu
217 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
218 fInnerRadius(rhs.fInnerRadius), fOuterRadi
219 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfL
220 fInnerStereo(rhs.fInnerStereo), fOuterSter
221 fTanInnerStereo(rhs.fTanInnerStereo), fTan
222 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInn
223 fOuterRadius2(rhs.fOuterRadius2), fTanInne
224 fTanOuterStereo2(rhs.fTanOuterStereo2),
225 fLowerEndcap(nullptr), fUpperEndcap(nullpt
226 fInnerHype(nullptr), fOuterHype(nullptr),
227 fCubicVolume(rhs.fCubicVolume), fSurfaceAr
228 {
229 for (auto i=0; i<2; ++i)
230 {
231 fEndZ[i] = rhs.fEndZ[i];
232 fEndInnerRadius[i] = rhs.fEndInnerRadius[i
233 fEndOuterRadius[i] = rhs.fEndOuterRadius[i
234 fEndPhi[i] = rhs.fEndPhi[i];
235 fEndZ2[i] = rhs.fEndZ2[i];
236 }
237 CreateSurfaces();
238 }
239
240
241 //============================================
242 //* Assignment operator ----------------------
243
244 G4TwistedTubs& G4TwistedTubs::operator = (cons
245 {
246 // Check assignment to self
247 //
248 if (this == &rhs) { return *this; }
249
250 // Copy base class data
251 //
252 G4VSolid::operator=(rhs);
253
254 // Copy data
255 //
256 fPhiTwist= rhs.fPhiTwist;
257 fInnerRadius= rhs.fInnerRadius; fOuterRadiu
258 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfL
259 fInnerStereo= rhs.fInnerStereo; fOuterStere
260 fTanInnerStereo= rhs.fTanInnerStereo; fTanO
261 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInn
262 fOuterRadius2= rhs.fOuterRadius2; fTanInner
263 fTanOuterStereo2= rhs.fTanOuterStereo2;
264 fLowerEndcap= fUpperEndcap= fLatterTwisted=
265 fInnerHype= fOuterHype= nullptr;
266 fCubicVolume= rhs.fCubicVolume; fSurfaceAre
267
268 for (auto i=0; i<2; ++i)
269 {
270 fEndZ[i] = rhs.fEndZ[i];
271 fEndInnerRadius[i] = rhs.fEndInnerRadius[
272 fEndOuterRadius[i] = rhs.fEndOuterRadius[
273 fEndPhi[i] = rhs.fEndPhi[i];
274 fEndZ2[i] = rhs.fEndZ2[i];
275 }
276
277 CreateSurfaces();
278 fRebuildPolyhedron = false;
279 delete fpPolyhedron; fpPolyhedron = nullptr
280
281 return *this;
282 }
283
284 //============================================
285 //* ComputeDimensions ------------------------
286
287 void G4TwistedTubs::ComputeDimensions(G4VPVPar
288 const G4
289 const G4
290 {
291 G4Exception("G4TwistedTubs::ComputeDimension
292 "GeomSolids0001", FatalException
293 "G4TwistedTubs does not support
294 }
295
296 //============================================
297 //* BoundingLimits ---------------------------
298
299 void G4TwistedTubs::BoundingLimits(G4ThreeVect
300 G4ThreeVect
301 {
302 // Find bounding tube
303 G4double rmin = GetInnerRadius();
304 G4double rmax = GetEndOuterRadius();
305
306 G4double zmin = std::min(GetEndZ(0), GetEndZ
307 G4double zmax = std::max(GetEndZ(0), GetEndZ
308
309 G4double dphi = 0.5*GetDPhi();
310 G4double sphi = std::min(GetEndPhi(0), GetEn
311 G4double ephi = std::max(GetEndPhi(0), GetEn
312 G4double totalphi = ephi - sphi;
313
314 // Find bounding box
315 if (dphi <= 0 || totalphi >= CLHEP::twopi)
316 {
317 pMin.set(-rmax,-rmax, zmin);
318 pMax.set( rmax, rmax, zmax);
319 }
320 else
321 {
322 G4TwoVector vmin,vmax;
323 G4GeomTools::DiskExtent(rmin, rmax, sphi,
324 pMin.set(vmin.x(), vmin.y(), zmin);
325 pMax.set(vmax.x(), vmax.y(), zmax);
326 }
327
328 // Check correctness of the bounding box
329 //
330 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
331 {
332 std::ostringstream message;
333 message << "Bad bounding box (min >= max)
334 << GetName() << " !"
335 << "\npMin = " << pMin
336 << "\npMax = " << pMax;
337 G4Exception("G4TwistedTubs::BoundingLimits
338 JustWarning, message);
339 DumpInfo();
340 }
341 }
342
343 //============================================
344 //* CalculateExtent --------------------------
345
346 G4bool
347 G4TwistedTubs::CalculateExtent(const EAxis pAx
348 const G4VoxelLi
349 const G4AffineT
350 G4double&
351 {
352 G4ThreeVector bmin, bmax;
353
354 // Get bounding box
355 BoundingLimits(bmin,bmax);
356
357 // Find extent
358 G4BoundingEnvelope bbox(bmin,bmax);
359 return bbox.CalculateExtent(pAxis,pVoxelLimi
360 }
361
362
363 //============================================
364 //* Inside -----------------------------------
365
366 EInside G4TwistedTubs::Inside(const G4ThreeVec
367 {
368
369 const G4double halftol
370 = 0.5 * G4GeometryTolerance::GetInstance(
371 // static G4int timerid = -1;
372 // G4Timer timer(timerid, "G4TwistedTubs",
373 // timer.Start();
374
375
376 EInside outerhypearea = ((G4TwistTubsHypeS
377 G4double innerhyperho = ((G4TwistTubsHypeS
378 G4double distanceToOut = p.getRho() - inner
379 EInside tmpinside;
380 if ((outerhypearea == kOutside) || (distanc
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 tmpinside;
401 }
402
403 //============================================
404 //* SurfaceNormal ----------------------------
405
406 G4ThreeVector G4TwistedTubs::SurfaceNormal(con
407 {
408 //
409 // return the normal unit vector to the Hyp
410 // p on (or nearly on) the surface
411 //
412 // Which of the three or four surfaces are
413 //
414
415
416 G4double distance = kInfinity;
417
418 G4VTwistSurface *surfaces[6];
419 surfaces[0] = fLatterTwisted;
420 surfaces[1] = fFormerTwisted;
421 surfaces[2] = fInnerHype;
422 surfaces[3] = fOuterHype;
423 surfaces[4] = fLowerEndcap;
424 surfaces[5] = fUpperEndcap;
425
426 G4ThreeVector xx;
427 G4ThreeVector bestxx;
428 G4int besti = -1;
429 for (auto i=0; i<6; ++i)
430 {
431 G4double tmpdistance = surfaces[i]->Dist
432 if (tmpdistance < distance)
433 {
434 distance = tmpdistance;
435 bestxx = xx;
436 besti = i;
437 }
438 }
439
440 return surfaces[besti]->GetNormal(bestxx, tr
441 }
442
443 //============================================
444 //* DistanceToIn (p, v) ----------------------
445
446 G4double G4TwistedTubs::DistanceToIn (const G4
447 const G4
448 {
449
450 // DistanceToIn (p, v):
451 // Calculate distance to surface of shape f
452 // along with the v, allowing for tolerance
453 // The function returns kInfinity if no int
454 // just grazing within tolerance.
455
456 //
457 // Calculate DistanceToIn(p,v)
458 //
459
460 EInside currentside = Inside(p);
461
462 if (currentside == kInside)
463 {
464 }
465 else
466 {
467 if (currentside == kSurface)
468 {
469 // particle is just on a boundary.
470 // If the particle is entering to the v
471 //
472 G4ThreeVector normal = SurfaceNormal(p)
473 if (normal*v < 0)
474 {
475 return 0;
476 }
477 }
478 }
479
480 // now, we can take smallest positive dista
481
482 // Initialize
483 //
484 G4double distance = kInfinity;
485
486 // find intersections and choose nearest on
487 //
488 G4VTwistSurface* surfaces[6];
489 surfaces[0] = fLowerEndcap;
490 surfaces[1] = fUpperEndcap;
491 surfaces[2] = fLatterTwisted;
492 surfaces[3] = fFormerTwisted;
493 surfaces[4] = fInnerHype;
494 surfaces[5] = fOuterHype;
495
496 G4ThreeVector xx;
497 G4ThreeVector bestxx;
498 for (const auto & surface : surfaces)
499 {
500 G4double tmpdistance = surface->Distance
501 if (tmpdistance < distance)
502 {
503 distance = tmpdistance;
504 bestxx = xx;
505 }
506 }
507 return distance;
508 }
509
510 //============================================
511 //* DistanceToIn (p) -------------------------
512
513 G4double G4TwistedTubs::DistanceToIn (const G4
514 {
515 // DistanceToIn(p):
516 // Calculate distance to surface of shape f
517 // allowing for tolerance
518
519 //
520 // Calculate DistanceToIn(p)
521 //
522
523 EInside currentside = Inside(p);
524
525 switch (currentside)
526 {
527 case (kInside) :
528 {}
529 case (kSurface) :
530 {
531 return 0;
532 }
533 case (kOutside) :
534 {
535 // Initialize
536 G4double distance = kInfinity;
537
538 // find intersections and choose near
539 G4VTwistSurface *surfaces[6];
540 surfaces[0] = fLowerEndcap;
541 surfaces[1] = fUpperEndcap;
542 surfaces[2] = fLatterTwisted;
543 surfaces[3] = fFormerTwisted;
544 surfaces[4] = fInnerHype;
545 surfaces[5] = fOuterHype;
546
547 G4ThreeVector xx;
548 G4ThreeVector bestxx;
549 for (const auto & surface : surfaces)
550 {
551 G4double tmpdistance = surface->Di
552 if (tmpdistance < distance)
553 {
554 distance = tmpdistance;
555 bestxx = xx;
556 }
557 }
558 return distance;
559 }
560 default :
561 {
562 G4Exception("G4TwistedTubs::DistanceT
563 FatalException, "Unknown
564 }
565 } // switch end
566
567 return kInfinity;
568 }
569
570 //============================================
571 //* DistanceToOut (p, v) ---------------------
572
573 G4double G4TwistedTubs::DistanceToOut( const G
574 const G
575 const G
576 G4bool
577 G4Three
578 {
579 // DistanceToOut (p, v):
580 // Calculate distance to surface of shape f
581 // along with the v, allowing for tolerance
582 // The function returns kInfinity if no int
583 // just grazing within tolerance.
584
585 //
586 // Calculate DistanceToOut(p,v)
587 //
588
589 EInside currentside = Inside(p);
590 if (currentside == kOutside)
591 {
592 }
593 else
594 {
595 if (currentside == kSurface)
596 {
597 // particle is just on a boundary.
598 // If the particle is exiting from the
599 //
600 G4ThreeVector normal = SurfaceNormal(p)
601 if (normal*v > 0)
602 {
603 if (calcNorm)
604 {
605 *norm = normal;
606 *validNorm = true;
607 }
608 return 0;
609 }
610 }
611 }
612
613 // now, we can take smallest positive dista
614
615 // Initialize
616 //
617 G4double distance = kInfinity;
618
619 // find intersections and choose nearest on
620 //
621 G4VTwistSurface* surfaces[6];
622 surfaces[0] = fLatterTwisted;
623 surfaces[1] = fFormerTwisted;
624 surfaces[2] = fInnerHype;
625 surfaces[3] = fOuterHype;
626 surfaces[4] = fLowerEndcap;
627 surfaces[5] = fUpperEndcap;
628
629 G4int besti = -1;
630 G4ThreeVector xx;
631 G4ThreeVector bestxx;
632 for (auto i=0; i<6; ++i)
633 {
634 G4double tmpdistance = surfaces[i]->Dist
635 if (tmpdistance < distance)
636 {
637 distance = tmpdistance;
638 bestxx = xx;
639 besti = i;
640 }
641 }
642
643 if (calcNorm)
644 {
645 if (besti != -1)
646 {
647 *norm = (surfaces[besti]->GetNormal(p
648 *validNorm = surfaces[besti]->IsValid
649 }
650 }
651
652 return distance;
653 }
654
655
656 //============================================
657 //* DistanceToOut (p) ------------------------
658
659 G4double G4TwistedTubs::DistanceToOut( const G
660 {
661 // DistanceToOut(p):
662 // Calculate distance to surface of shape f
663 // allowing for tolerance
664
665 //
666 // Calculate DistanceToOut(p)
667 //
668
669 EInside currentside = Inside(p);
670
671 switch (currentside)
672 {
673 case (kOutside) :
674 {
675 }
676 case (kSurface) :
677 {
678 return 0;
679 }
680 case (kInside) :
681 {
682 // Initialize
683 G4double distance = kInfinity;
684
685 // find intersections and choose near
686 G4VTwistSurface* surfaces[6];
687 surfaces[0] = fLatterTwisted;
688 surfaces[1] = fFormerTwisted;
689 surfaces[2] = fInnerHype;
690 surfaces[3] = fOuterHype;
691 surfaces[4] = fLowerEndcap;
692 surfaces[5] = fUpperEndcap;
693
694 G4ThreeVector xx;
695 G4ThreeVector bestxx;
696 for (const auto & surface : surfaces)
697 {
698 G4double tmpdistance = surface->Di
699 if (tmpdistance < distance)
700 {
701 distance = tmpdistance;
702 bestxx = xx;
703 }
704 }
705 return distance;
706 }
707 default :
708 {
709 G4Exception("G4TwistedTubs::DistanceT
710 FatalException, "Unknown
711 }
712 } // switch end
713
714 return 0.;
715 }
716
717 //============================================
718 //* StreamInfo -------------------------------
719
720 std::ostream& G4TwistedTubs::StreamInfo(std::o
721 {
722 //
723 // Stream object contents to an output strea
724 //
725 G4long oldprc = os.precision(16);
726 os << "-------------------------------------
727 << " *** Dump for solid - " << GetName
728 << " =================================
729 << " Solid type: G4TwistedTubs\n"
730 << " Parameters: \n"
731 << " -ve end Z : " << fEn
732 << " +ve end Z : " << fEn
733 << " inner end radius(-ve z): " << fEn
734 << " inner end radius(+ve z): " << fEn
735 << " outer end radius(-ve z): " << fEn
736 << " outer end radius(+ve z): " << fEn
737 << " inner radius (z=0) : " << fIn
738 << " outer radius (z=0) : " << fOu
739 << " twisted angle : " << fPh
740 << " inner stereo angle : " << fIn
741 << " outer stereo angle : " << fOu
742 << " phi-width of a piece : " << fDP
743 << "-------------------------------------
744 os.precision(oldprc);
745
746 return os;
747 }
748
749
750 //============================================
751 //* DiscribeYourselfTo -----------------------
752
753 void G4TwistedTubs::DescribeYourselfTo (G4VGra
754 {
755 scene.AddSolid (*this);
756 }
757
758 //============================================
759 //* GetExtent --------------------------------
760
761 G4VisExtent G4TwistedTubs::GetExtent() const
762 {
763 // Define the sides of the box into which th
764 //
765 G4ThreeVector pmin,pmax;
766 BoundingLimits(pmin,pmax);
767 return { pmin.x(),pmax.x(),
768 pmin.y(),pmax.y(),
769 pmin.z(),pmax.z() };
770 }
771
772 //============================================
773 //* CreatePolyhedron -------------------------
774
775 G4Polyhedron* G4TwistedTubs::CreatePolyhedron
776 {
777 // number of meshes
778 //
779 G4double absPhiTwist = std::abs(fPhiTwist);
780 G4double dA = std::max(fDPhi,absPhiTwist);
781 const G4int k =
782 G4int(G4Polyhedron::GetNumberOfRotationSte
783 const G4int n =
784 G4int(G4Polyhedron::GetNumberOfRotationSte
785
786 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
787 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)
788
789 auto ph = new G4Polyhedron;
790 typedef G4double G4double3[3];
791 typedef G4int G4int4[4];
792 auto xyz = new G4double3[nnodes]; // number
793 auto faces = new G4int4[nfaces] ; // number
794 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
795 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
796 fInnerHype->GetFacets(k,n,xyz,faces,2) ;
797 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
798 fOuterHype->GetFacets(k,n,xyz,faces,4) ;
799 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
800
801 ph->createPolyhedron(nnodes,nfaces,xyz,faces
802
803 delete[] xyz;
804 delete[] faces;
805
806 return ph;
807 }
808
809 //============================================
810 //* GetPolyhedron ----------------------------
811
812 G4Polyhedron* G4TwistedTubs::GetPolyhedron ()
813 {
814 if (fpPolyhedron == nullptr ||
815 fRebuildPolyhedron ||
816 fpPolyhedron->GetNumberOfRotationStepsAt
817 fpPolyhedron->GetNumberOfRotationSteps()
818 {
819 G4AutoLock l(&polyhedronMutex);
820 delete fpPolyhedron;
821 fpPolyhedron = CreatePolyhedron();
822 fRebuildPolyhedron = false;
823 l.unlock();
824 }
825 return fpPolyhedron;
826 }
827
828 //============================================
829 //* CreateSurfaces ---------------------------
830
831 void G4TwistedTubs::CreateSurfaces()
832 {
833 // create 6 surfaces of TwistedTub
834
835 G4ThreeVector x0(0, 0, fEndZ[0]);
836 G4ThreeVector n (0, 0, -1);
837
838 fLowerEndcap = new G4TwistTubsFlatSide("Low
839 fEndInnerR
840 fDPhi, fEn
841
842 fUpperEndcap = new G4TwistTubsFlatSide("Upp
843 fEndInnerR
844 fDPhi, fEn
845
846 G4RotationMatrix rotHalfDPhi;
847 rotHalfDPhi.rotateZ(0.5*fDPhi);
848
849 fLatterTwisted = new G4TwistTubsSide("Latte
850 fEndI
851 fDPhi
852 fInne
853 1 ) ;
854 fFormerTwisted = new G4TwistTubsSide("Forme
855 fEndI
856 fDPhi
857 fInne
858 -1 )
859
860 fInnerHype = new G4TwistTubsHypeSide("Inner
861 fEndIn
862 fDPhi,
863 fInner
864 fTanIn
865 fOuterHype = new G4TwistTubsHypeSide("Outer
866 fEndIn
867 fDPhi,
868 fInner
869 fTanIn
870
871
872 // set neighbour surfaces
873 //
874 fLowerEndcap->SetNeighbours(fInnerHype, fLa
875 fOuterHype, fFo
876 fUpperEndcap->SetNeighbours(fInnerHype, fLa
877 fOuterHype, fFo
878 fLatterTwisted->SetNeighbours(fInnerHype, f
879 fOuterHype, f
880 fFormerTwisted->SetNeighbours(fInnerHype, f
881 fOuterHype, f
882 fInnerHype->SetNeighbours(fLatterTwisted, f
883 fFormerTwisted, f
884 fOuterHype->SetNeighbours(fLatterTwisted, f
885 fFormerTwisted, f
886 }
887
888
889 //============================================
890 //* GetEntityType ----------------------------
891
892 G4GeometryType G4TwistedTubs::GetEntityType()
893 {
894 return {"G4TwistedTubs"};
895 }
896
897 //============================================
898 //* Clone ------------------------------------
899
900 G4VSolid* G4TwistedTubs::Clone() const
901 {
902 return new G4TwistedTubs(*this);
903 }
904
905 //============================================
906 //* GetCubicVolume ---------------------------
907
908 G4double G4TwistedTubs::GetCubicVolume()
909 {
910 if (fCubicVolume == 0.)
911 {
912 G4double DPhi = GetDPhi();
913 G4double Z0 = GetEndZ(0);
914 G4double Z1 = GetEndZ(1);
915 G4double Ain = GetInnerRadius();
916 G4double Aout = GetOuterRadius();
917 G4double R0in = GetEndInnerRadius(0);
918 G4double R1in = GetEndInnerRadius(1);
919 G4double R0out = GetEndOuterRadius(0);
920 G4double R1out = GetEndOuterRadius(1);
921
922 // V_hyperboloid = pi*h*(2*a*a + R*R)/3
923 fCubicVolume = (2.*(Z1 - Z0)*(Aout + Ain)*
924 + Z1*(R1out + R1in)*(R1out
925 - Z0*(R0out + R0in)*(R0out
926 }
927 return fCubicVolume;
928 }
929
930 //============================================
931 //* GetLateralArea ---------------------------
932
933 G4double
934 G4TwistedTubs::GetLateralArea(G4double a, G4do
935 {
936 if (z == 0) return 0.;
937 G4double h = std::abs(z);
938 G4double area = h*a;
939 if (std::abs(a - r) > kCarTolerance)
940 {
941 G4double aa = a*a;
942 G4double hh = h*h;
943 G4double rr = r*r;
944 G4double cc = aa*hh/(rr - aa);
945 G4double k = std::sqrt(aa + cc)/cc;
946 G4double kh = k*h;
947 area = 0.5*a*(h*std::sqrt(1. + kh*kh) + st
948 }
949 return GetDPhi()*area;
950 }
951
952 //============================================
953 //* GetPhiCutArea ----------------------------
954
955 G4double
956 G4TwistedTubs::GetPhiCutArea(G4double a, G4dou
957 {
958 if (GetDPhi() >= CLHEP::twopi || r <= 0 || z
959 G4double h = std::abs(z);
960 G4double area = h*a;
961 if (GetPhiTwist() > kCarTolerance)
962 {
963 G4double sinw = std::sin(0.5*GetPhiTwist()
964 G4double p = sinw*r/h;
965 G4double q = sinw*r/a;
966 G4double pp = p*p;
967 G4double qq = q*q;
968 G4double pq = p*q;
969 G4double sqroot = std::sqrt(pp + qq + 1);
970 area = (pq*sqroot +
971 0.5*p*(pp + 3.)*std::atanh(q/sqroo
972 0.5*q*(qq + 3.)*std::atanh(p/sqroo
973 std::atan(sqroot/(pq)) - CLHEP::ha
974 }
975 return area;
976 }
977
978 //============================================
979 //* GetSurfaceArea ---------------------------
980
981 G4double G4TwistedTubs::GetSurfaceArea()
982 {
983 if (fSurfaceArea == 0.)
984 {
985 G4double dphi = GetDPhi();
986 G4double Ainn = GetInnerRadius();
987 G4double Aout = GetOuterRadius();
988 G4double Rinn0 = GetEndInnerRadius(0);
989 G4double Rout0 = GetEndOuterRadius(0);
990 G4double Rinn1 = GetEndInnerRadius(1);
991 G4double Rout1 = GetEndOuterRadius(1);
992 G4double z0 = GetEndZ(0);
993 G4double z1 = GetEndZ(1);
994
995 G4double base0 = 0.5*dphi*(Rout0*Rout0 - R
996 G4double inner0 = GetLateralArea(Ainn, Rin
997 G4double outer0 = GetLateralArea(Aout, Rou
998 G4double cut0 =
999 GetPhiCutArea(Aout, Rout0, z0) - GetPhiC
1000
1001 G4double base1 = base0;
1002 G4double inner1 = inner0;
1003 G4double outer1 = outer0;
1004 G4double cut1 = cut0;
1005 if (std::abs(z0) != std::abs(z1))
1006 {
1007 base1 = 0.5*dphi*(Rout1*Rout1 - Rinn1*R
1008 inner1 = GetLateralArea(Ainn, Rinn1, z1
1009 outer1 = GetLateralArea(Aout, Rout1, z1
1010 cut1 =
1011 GetPhiCutArea(Aout, Rout1, z1) - GetPhi
1012 }
1013 fSurfaceArea = base0 + base1 +
1014 ((z0*z1 < 0) ?
1015 (inner0 + inner1 + outer0 + outer1 + 2.
1016 std::abs(inner0 - inner1 + outer0 - out
1017 }
1018 return fSurfaceArea;
1019 }
1020
1021 //===========================================
1022 //* GetPointOnSurface -----------------------
1023
1024 G4ThreeVector G4TwistedTubs::GetPointOnSurfac
1025 {
1026
1027 G4double z = G4RandFlat::shoot(fEndZ[0],fEn
1028 G4double phi , phimin, phimax ;
1029 G4double x , xmin, xmax ;
1030 G4double r , rmin, rmax ;
1031
1032 G4double a1 = fOuterHype->GetSurfaceArea()
1033 G4double a2 = fInnerHype->GetSurfaceArea()
1034 G4double a3 = fLatterTwisted->GetSurfaceAre
1035 G4double a4 = fFormerTwisted->GetSurfaceAre
1036 G4double a5 = fLowerEndcap->GetSurfaceArea(
1037 G4double a6 = fUpperEndcap->GetSurfaceArea(
1038
1039 G4double chose = G4RandFlat::shoot(0.,a1 +
1040
1041 if(chose < a1)
1042 {
1043
1044 phimin = fOuterHype->GetBoundaryMin(z) ;
1045 phimax = fOuterHype->GetBoundaryMax(z) ;
1046 phi = G4RandFlat::shoot(phimin,phimax) ;
1047
1048 return fOuterHype->SurfacePoint(phi,z,tru
1049
1050 }
1051 else if ( (chose >= a1) && (chose < a1 + a2
1052 {
1053
1054 phimin = fInnerHype->GetBoundaryMin(z) ;
1055 phimax = fInnerHype->GetBoundaryMax(z) ;
1056 phi = G4RandFlat::shoot(phimin,phimax) ;
1057
1058 return fInnerHype->SurfacePoint(phi,z,tru
1059
1060 }
1061 else if ( (chose >= a1 + a2 ) && (chose < a
1062 {
1063
1064 xmin = fLatterTwisted->GetBoundaryMin(z)
1065 xmax = fLatterTwisted->GetBoundaryMax(z)
1066 x = G4RandFlat::shoot(xmin,xmax) ;
1067
1068 return fLatterTwisted->SurfacePoint(x,z,t
1069
1070 }
1071 else if ( (chose >= a1 + a2 + a3 ) && (cho
1072 {
1073
1074 xmin = fFormerTwisted->GetBoundaryMin(z)
1075 xmax = fFormerTwisted->GetBoundaryMax(z)
1076 x = G4RandFlat::shoot(xmin,xmax) ;
1077
1078 return fFormerTwisted->SurfacePoint(x,z,t
1079 }
1080 else if( (chose >= a1 + a2 + a3 + a4 )&&(c
1081 {
1082 rmin = GetEndInnerRadius(0) ;
1083 rmax = GetEndOuterRadius(0) ;
1084 r = std::sqrt(G4RandFlat::shoot()*(sqr(rm
1085
1086 phimin = fLowerEndcap->GetBoundaryMin(r)
1087 phimax = fLowerEndcap->GetBoundaryMax(r)
1088 phi = G4RandFlat::shoot(phimin,phimax)
1089
1090 return fLowerEndcap->SurfacePoint(phi,r,t
1091 }
1092 else
1093 {
1094 rmin = GetEndInnerRadius(1) ;
1095 rmax = GetEndOuterRadius(1) ;
1096 r = rmin + (rmax-rmin)*std::sqrt(G4RandFl
1097
1098 phimin = fUpperEndcap->GetBoundaryMin(r)
1099 phimax = fUpperEndcap->GetBoundaryMax(r)
1100 phi = G4RandFlat::shoot(phimin,phimax)
1101
1102 return fUpperEndcap->SurfacePoint(phi,r,t
1103 }
1104 }
1105