| Geant4 Cross Reference |
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 (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupiter-2.5.02 application.
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_INITIALIZER;
57 }
58
59
60 //=====================================================================
61 //* constructors ------------------------------------------------------
62
63 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
64 G4double twistedangle,
65 G4double endinnerrad,
66 G4double endouterrad,
67 G4double halfzlen,
68 G4double dphi)
69 : G4VSolid(pname), fDPhi(dphi),
70 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
71 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
72 {
73 if (endinnerrad < DBL_MIN)
74 {
75 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
76 FatalErrorInArgument, "Invalid end-inner-radius!");
77 }
78
79 G4double sinhalftwist = std::sin(0.5 * twistedangle);
80
81 G4double endinnerradX = endinnerrad * sinhalftwist;
82 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
83 - endinnerradX * endinnerradX );
84
85 G4double endouterradX = endouterrad * sinhalftwist;
86 G4double outerrad = std::sqrt( endouterrad * endouterrad
87 - endouterradX * endouterradX );
88
89 // temporary treatment!!
90 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
91 CreateSurfaces();
92 }
93
94 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
95 G4double twistedangle,
96 G4double endinnerrad,
97 G4double endouterrad,
98 G4double halfzlen,
99 G4int nseg,
100 G4double totphi)
101 : G4VSolid(pname),
102 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
103 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
104 {
105
106 if (nseg == 0)
107 {
108 std::ostringstream message;
109 message << "Invalid number of segments." << G4endl
110 << " nseg = " << nseg;
111 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
112 FatalErrorInArgument, message);
113 }
114 if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
115 {
116 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
117 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
118 }
119
120 G4double sinhalftwist = std::sin(0.5 * twistedangle);
121
122 G4double endinnerradX = endinnerrad * sinhalftwist;
123 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
124 - endinnerradX * endinnerradX );
125
126 G4double endouterradX = endouterrad * sinhalftwist;
127 G4double outerrad = std::sqrt( endouterrad * endouterrad
128 - endouterradX * endouterradX );
129
130 // temporary treatment!!
131 fDPhi = totphi / nseg;
132 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
133 CreateSurfaces();
134 }
135
136 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
137 G4double twistedangle,
138 G4double innerrad,
139 G4double outerrad,
140 G4double negativeEndz,
141 G4double positiveEndz,
142 G4double dphi)
143 : G4VSolid(pname), fDPhi(dphi),
144 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
145 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
146 {
147 if (innerrad < DBL_MIN)
148 {
149 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
150 FatalErrorInArgument, "Invalid end-inner-radius!");
151 }
152
153 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
154 CreateSurfaces();
155 }
156
157 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
158 G4double twistedangle,
159 G4double innerrad,
160 G4double outerrad,
161 G4double negativeEndz,
162 G4double positiveEndz,
163 G4int nseg,
164 G4double totphi)
165 : G4VSolid(pname),
166 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
167 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
168 {
169 if (nseg == 0)
170 {
171 std::ostringstream message;
172 message << "Invalid number of segments." << G4endl
173 << " nseg = " << nseg;
174 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
175 FatalErrorInArgument, message);
176 }
177 if (totphi == DBL_MIN || innerrad < DBL_MIN)
178 {
179 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
180 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
181 }
182
183 fDPhi = totphi / nseg;
184 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
185 CreateSurfaces();
186 }
187
188 //=====================================================================
189 //* Fake default constructor ------------------------------------------
190
191 G4TwistedTubs::G4TwistedTubs( __void__& a )
192 : G4VSolid(a),
193 fLowerEndcap(nullptr), fUpperEndcap(nullptr),
194 fLatterTwisted(nullptr), fFormerTwisted(nullptr),
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 G4TwistedTubs& rhs)
217 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
218 fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
219 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
220 fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
221 fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
222 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
223 fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
224 fTanOuterStereo2(rhs.fTanOuterStereo2),
225 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr), fFormerTwisted(nullptr),
226 fInnerHype(nullptr), fOuterHype(nullptr),
227 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea)
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 = (const G4TwistedTubs& rhs)
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; fOuterRadius= rhs.fOuterRadius;
258 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
259 fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
260 fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
261 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
262 fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
263 fTanOuterStereo2= rhs.fTanOuterStereo2;
264 fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= nullptr;
265 fInnerHype= fOuterHype= nullptr;
266 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
267
268 for (auto i=0; i<2; ++i)
269 {
270 fEndZ[i] = rhs.fEndZ[i];
271 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
272 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
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(G4VPVParameterisation* /* p */ ,
288 const G4int /* n */ ,
289 const G4VPhysicalVolume* /* pRep */ )
290 {
291 G4Exception("G4TwistedTubs::ComputeDimensions()",
292 "GeomSolids0001", FatalException,
293 "G4TwistedTubs does not support Parameterisation.");
294 }
295
296 //=====================================================================
297 //* BoundingLimits ----------------------------------------------------
298
299 void G4TwistedTubs::BoundingLimits(G4ThreeVector& pMin,
300 G4ThreeVector& pMax) const
301 {
302 // Find bounding tube
303 G4double rmin = GetInnerRadius();
304 G4double rmax = GetEndOuterRadius();
305
306 G4double zmin = std::min(GetEndZ(0), GetEndZ(1));
307 G4double zmax = std::max(GetEndZ(0), GetEndZ(1));
308
309 G4double dphi = 0.5*GetDPhi();
310 G4double sphi = std::min(GetEndPhi(0), GetEndPhi(1)) - dphi;
311 G4double ephi = std::max(GetEndPhi(0), GetEndPhi(1)) + dphi;
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, totalphi, vmin, vmax);
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.y() || pMin.z() >= pMax.z())
331 {
332 std::ostringstream message;
333 message << "Bad bounding box (min >= max) for solid: "
334 << GetName() << " !"
335 << "\npMin = " << pMin
336 << "\npMax = " << pMax;
337 G4Exception("G4TwistedTubs::BoundingLimits()", "GeomMgt0001",
338 JustWarning, message);
339 DumpInfo();
340 }
341 }
342
343 //=====================================================================
344 //* CalculateExtent ---------------------------------------------------
345
346 G4bool
347 G4TwistedTubs::CalculateExtent(const EAxis pAxis,
348 const G4VoxelLimits& pVoxelLimit,
349 const G4AffineTransform& pTransform,
350 G4double& pMin, G4double& pMax) const
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,pVoxelLimit,pTransform,pMin,pMax);
360 }
361
362
363 //=====================================================================
364 //* Inside ------------------------------------------------------------
365
366 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
367 {
368
369 const G4double halftol
370 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
371 // static G4int timerid = -1;
372 // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
373 // timer.Start();
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 EInside tmpinside;
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 tmpinside;
401 }
402
403 //=====================================================================
404 //* SurfaceNormal -----------------------------------------------------
405
406 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
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
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]->DistanceTo(p, xx);
432 if (tmpdistance < distance)
433 {
434 distance = tmpdistance;
435 bestxx = xx;
436 besti = i;
437 }
438 }
439
440 return surfaces[besti]->GetNormal(bestxx, true);
441 }
442
443 //=====================================================================
444 //* DistanceToIn (p, v) -----------------------------------------------
445
446 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
447 const G4ThreeVector& v ) const
448 {
449
450 // DistanceToIn (p, v):
451 // Calculate distance to surface of shape from `outside'
452 // along with the v, allowing for tolerance.
453 // The function returns kInfinity if no intersection or
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 volume, return 0.
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 distance.
481
482 // Initialize
483 //
484 G4double distance = kInfinity;
485
486 // find intersections and choose nearest one.
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->DistanceToIn(p, v, xx);
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 G4ThreeVector& p) const
514 {
515 // DistanceToIn(p):
516 // Calculate distance to surface of shape from `outside',
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 nearest one.
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->DistanceTo(p, xx);
552 if (tmpdistance < distance)
553 {
554 distance = tmpdistance;
555 bestxx = xx;
556 }
557 }
558 return distance;
559 }
560 default :
561 {
562 G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
563 FatalException, "Unknown point location!");
564 }
565 } // switch end
566
567 return kInfinity;
568 }
569
570 //=====================================================================
571 //* DistanceToOut (p, v) ----------------------------------------------
572
573 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
574 const G4ThreeVector& v,
575 const G4bool calcNorm,
576 G4bool *validNorm,
577 G4ThreeVector *norm ) const
578 {
579 // DistanceToOut (p, v):
580 // Calculate distance to surface of shape from `inside'
581 // along with the v, allowing for tolerance.
582 // The function returns kInfinity if no intersection or
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 volume, return 0.
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 distance.
614
615 // Initialize
616 //
617 G4double distance = kInfinity;
618
619 // find intersections and choose nearest one.
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]->DistanceToOut(p, v, xx);
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, true));
648 *validNorm = surfaces[besti]->IsValidNorm();
649 }
650 }
651
652 return distance;
653 }
654
655
656 //=====================================================================
657 //* DistanceToOut (p) ----------------------------------------------
658
659 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
660 {
661 // DistanceToOut(p):
662 // Calculate distance to surface of shape from `inside',
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 nearest one.
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->DistanceTo(p, xx);
699 if (tmpdistance < distance)
700 {
701 distance = tmpdistance;
702 bestxx = xx;
703 }
704 }
705 return distance;
706 }
707 default :
708 {
709 G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
710 FatalException, "Unknown point location!");
711 }
712 } // switch end
713
714 return 0.;
715 }
716
717 //=====================================================================
718 //* StreamInfo --------------------------------------------------------
719
720 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
721 {
722 //
723 // Stream object contents to an output stream
724 //
725 G4long oldprc = os.precision(16);
726 os << "-----------------------------------------------------------\n"
727 << " *** Dump for solid - " << GetName() << " ***\n"
728 << " ===================================================\n"
729 << " Solid type: G4TwistedTubs\n"
730 << " Parameters: \n"
731 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
732 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
733 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
734 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
735 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
736 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
737 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
738 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
739 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
740 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
741 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
742 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
743 << "-----------------------------------------------------------\n";
744 os.precision(oldprc);
745
746 return os;
747 }
748
749
750 //=====================================================================
751 //* DiscribeYourselfTo ------------------------------------------------
752
753 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
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 the G4Tubs instance would fit.
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 () const
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::GetNumberOfRotationSteps() * dA / twopi) + 2;
783 const G4int n =
784 G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
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)*(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 of nodes
793 auto faces = new G4int4[nfaces] ; // number of faces
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 () const
813 {
814 if (fpPolyhedron == nullptr ||
815 fRebuildPolyhedron ||
816 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
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("LowerEndcap",
839 fEndInnerRadius, fEndOuterRadius,
840 fDPhi, fEndPhi, fEndZ, -1) ;
841
842 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
843 fEndInnerRadius, fEndOuterRadius,
844 fDPhi, fEndPhi, fEndZ, 1) ;
845
846 G4RotationMatrix rotHalfDPhi;
847 rotHalfDPhi.rotateZ(0.5*fDPhi);
848
849 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
850 fEndInnerRadius, fEndOuterRadius,
851 fDPhi, fEndPhi, fEndZ,
852 fInnerRadius, fOuterRadius, fKappa,
853 1 ) ;
854 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
855 fEndInnerRadius, fEndOuterRadius,
856 fDPhi, fEndPhi, fEndZ,
857 fInnerRadius, fOuterRadius, fKappa,
858 -1 ) ;
859
860 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
861 fEndInnerRadius, fEndOuterRadius,
862 fDPhi, fEndPhi, fEndZ,
863 fInnerRadius, fOuterRadius,fKappa,
864 fTanInnerStereo, fTanOuterStereo, -1) ;
865 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
866 fEndInnerRadius, fEndOuterRadius,
867 fDPhi, fEndPhi, fEndZ,
868 fInnerRadius, fOuterRadius,fKappa,
869 fTanInnerStereo, fTanOuterStereo, 1) ;
870
871
872 // set neighbour surfaces
873 //
874 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
875 fOuterHype, fFormerTwisted);
876 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
877 fOuterHype, fFormerTwisted);
878 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
879 fOuterHype, fUpperEndcap);
880 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
881 fOuterHype, fUpperEndcap);
882 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
883 fFormerTwisted, fUpperEndcap);
884 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
885 fFormerTwisted, fUpperEndcap);
886 }
887
888
889 //=====================================================================
890 //* GetEntityType -----------------------------------------------------
891
892 G4GeometryType G4TwistedTubs::GetEntityType() const
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)*(Aout - Ain)
924 + Z1*(R1out + R1in)*(R1out - R1in)
925 - Z0*(R0out + R0in)*(R0out - R0in))*DPhi/6.;
926 }
927 return fCubicVolume;
928 }
929
930 //=====================================================================
931 //* GetLateralArea ----------------------------------------------------
932
933 G4double
934 G4TwistedTubs::GetLateralArea(G4double a, G4double r, G4double z) const
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) + std::asinh(kh)/k);
948 }
949 return GetDPhi()*area;
950 }
951
952 //=====================================================================
953 //* GetPhiCutArea -----------------------------------------------------
954
955 G4double
956 G4TwistedTubs::GetPhiCutArea(G4double a, G4double r, G4double z) const
957 {
958 if (GetDPhi() >= CLHEP::twopi || r <= 0 || z == 0) return 0.;
959 G4double h = std::abs(z);
960 G4double area = h*a;
961 if (GetPhiTwist() > kCarTolerance)
962 {
963 G4double sinw = std::sin(0.5*GetPhiTwist())*h/GetZHalfLength();
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/sqroot) +
972 0.5*q*(qq + 3.)*std::atanh(p/sqroot) +
973 std::atan(sqroot/(pq)) - CLHEP::halfpi)*h*a/(3.*pq);
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 - Rinn0*Rinn0); // lower base
996 G4double inner0 = GetLateralArea(Ainn, Rinn0, z0); // lower inner surface
997 G4double outer0 = GetLateralArea(Aout, Rout0, z0); // lower outer surface
998 G4double cut0 = // lower phi cut
999 GetPhiCutArea(Aout, Rout0, z0) - GetPhiCutArea(Ainn, Rinn0, z0);
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*Rinn1); // upper base
1008 inner1 = GetLateralArea(Ainn, Rinn1, z1); // upper inner surface
1009 outer1 = GetLateralArea(Aout, Rout1, z1); // upper outer surface
1010 cut1 = // upper phi cut
1011 GetPhiCutArea(Aout, Rout1, z1) - GetPhiCutArea(Ainn, Rinn1, z1);
1012 }
1013 fSurfaceArea = base0 + base1 +
1014 ((z0*z1 < 0) ?
1015 (inner0 + inner1 + outer0 + outer1 + 2.*(cut0 + cut1)) :
1016 std::abs(inner0 - inner1 + outer0 - outer1 + 2.*(cut0 - cut1)));
1017 }
1018 return fSurfaceArea;
1019 }
1020
1021 //=====================================================================
1022 //* GetPointOnSurface -------------------------------------------------
1023
1024 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1025 {
1026
1027 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
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->GetSurfaceArea() ;
1035 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1036 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1037 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1038
1039 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
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,true) ;
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,true) ;
1059
1060 }
1061 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
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,true) ;
1069
1070 }
1071 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
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,true) ;
1079 }
1080 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1081 {
1082 rmin = GetEndInnerRadius(0) ;
1083 rmax = GetEndOuterRadius(0) ;
1084 r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1085
1086 phimin = fLowerEndcap->GetBoundaryMin(r) ;
1087 phimax = fLowerEndcap->GetBoundaryMax(r) ;
1088 phi = G4RandFlat::shoot(phimin,phimax) ;
1089
1090 return fLowerEndcap->SurfacePoint(phi,r,true) ;
1091 }
1092 else
1093 {
1094 rmin = GetEndInnerRadius(1) ;
1095 rmax = GetEndOuterRadius(1) ;
1096 r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1097
1098 phimin = fUpperEndcap->GetBoundaryMin(r) ;
1099 phimax = fUpperEndcap->GetBoundaryMax(r) ;
1100 phi = G4RandFlat::shoot(phimin,phimax) ;
1101
1102 return fUpperEndcap->SurfacePoint(phi,r,true) ;
1103 }
1104 }
1105