Geant4 11.2.2
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
35#include "G4GeomTools.hh"
37#include "G4SystemOfUnits.hh"
39#include "G4VoxelLimits.hh"
40#include "G4AffineTransform.hh"
41#include "G4BoundingEnvelope.hh"
42#include "G4ClippablePolygon.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
54namespace
55{
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
57}
58
59//=====================================================================
60//* constructors ------------------------------------------------------
61
63 G4double twistedangle,
64 G4double endinnerrad,
65 G4double endouterrad,
66 G4double halfzlen,
67 G4double dphi)
68 : G4VSolid(pname), fDPhi(dphi),
69 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
70 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
71{
72 if (endinnerrad < DBL_MIN)
73 {
74 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
75 FatalErrorInArgument, "Invalid end-inner-radius!");
76 }
77
78 G4double sinhalftwist = std::sin(0.5 * twistedangle);
79
80 G4double endinnerradX = endinnerrad * sinhalftwist;
81 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
82 - endinnerradX * endinnerradX );
83
84 G4double endouterradX = endouterrad * sinhalftwist;
85 G4double outerrad = std::sqrt( endouterrad * endouterrad
86 - endouterradX * endouterradX );
87
88 // temporary treatment!!
89 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
90 CreateSurfaces();
91}
92
94 G4double twistedangle,
95 G4double endinnerrad,
96 G4double endouterrad,
97 G4double halfzlen,
98 G4int nseg,
99 G4double totphi)
100 : G4VSolid(pname),
101 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
102 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
103{
104
105 if (nseg == 0)
106 {
107 std::ostringstream message;
108 message << "Invalid number of segments." << G4endl
109 << " nseg = " << nseg;
110 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
111 FatalErrorInArgument, message);
112 }
113 if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
114 {
115 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
116 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
117 }
118
119 G4double sinhalftwist = std::sin(0.5 * twistedangle);
120
121 G4double endinnerradX = endinnerrad * sinhalftwist;
122 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
123 - endinnerradX * endinnerradX );
124
125 G4double endouterradX = endouterrad * sinhalftwist;
126 G4double outerrad = std::sqrt( endouterrad * endouterrad
127 - endouterradX * endouterradX );
128
129 // temporary treatment!!
130 fDPhi = totphi / nseg;
131 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
132 CreateSurfaces();
133}
134
136 G4double twistedangle,
137 G4double innerrad,
138 G4double outerrad,
139 G4double negativeEndz,
140 G4double positiveEndz,
141 G4double dphi)
142 : G4VSolid(pname), fDPhi(dphi),
143 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
144 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
145{
146 if (innerrad < DBL_MIN)
147 {
148 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
149 FatalErrorInArgument, "Invalid end-inner-radius!");
150 }
151
152 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
153 CreateSurfaces();
154}
155
157 G4double twistedangle,
158 G4double innerrad,
159 G4double outerrad,
160 G4double negativeEndz,
161 G4double positiveEndz,
162 G4int nseg,
163 G4double totphi)
164 : G4VSolid(pname),
165 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
166 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
167{
168 if (nseg == 0)
169 {
170 std::ostringstream message;
171 message << "Invalid number of segments." << G4endl
172 << " nseg = " << nseg;
173 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
174 FatalErrorInArgument, message);
175 }
176 if (totphi == DBL_MIN || innerrad < DBL_MIN)
177 {
178 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
179 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
180 }
181
182 fDPhi = totphi / nseg;
183 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
184 CreateSurfaces();
185}
186
187//=====================================================================
188//* Fake default constructor ------------------------------------------
189
191 : G4VSolid(a),
192 fLowerEndcap(nullptr), fUpperEndcap(nullptr),
193 fLatterTwisted(nullptr), fFormerTwisted(nullptr),
194 fInnerHype(nullptr), fOuterHype(nullptr)
195{
196}
197
198//=====================================================================
199//* destructor --------------------------------------------------------
200
202{
203 delete fLowerEndcap;
204 delete fUpperEndcap;
205 delete fLatterTwisted;
206 delete fFormerTwisted;
207 delete fInnerHype;
208 delete fOuterHype;
209 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(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr), fFormerTwisted(nullptr),
225 fInnerHype(nullptr), fOuterHype(nullptr),
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= nullptr;
269 fInnerHype= fOuterHype= nullptr;
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 // Find bounding tube
312 G4double rmin = GetInnerRadius();
314
315 G4double zmin = std::min(GetEndZ(0), GetEndZ(1));
316 G4double zmax = std::max(GetEndZ(0), GetEndZ(1));
317
318 G4double dphi = 0.5*GetDPhi();
319 G4double sphi = std::min(GetEndPhi(0), GetEndPhi(1)) - dphi;
320 G4double ephi = std::max(GetEndPhi(0), GetEndPhi(1)) + dphi;
321 G4double totalphi = ephi - sphi;
322
323 // Find bounding box
324 if (dphi <= 0 || totalphi >= CLHEP::twopi)
325 {
326 pMin.set(-rmax,-rmax, zmin);
327 pMax.set( rmax, rmax, zmax);
328 }
329 else
330 {
331 G4TwoVector vmin,vmax;
332 G4GeomTools::DiskExtent(rmin, rmax, sphi, totalphi, vmin, vmax);
333 pMin.set(vmin.x(), vmin.y(), zmin);
334 pMax.set(vmax.x(), vmax.y(), zmax);
335 }
336
337 // Check correctness of the bounding box
338 //
339 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
340 {
341 std::ostringstream message;
342 message << "Bad bounding box (min >= max) for solid: "
343 << GetName() << " !"
344 << "\npMin = " << pMin
345 << "\npMax = " << pMax;
346 G4Exception("G4TwistedTubs::BoundingLimits()", "GeomMgt0001",
347 JustWarning, message);
348 DumpInfo();
349 }
350}
351
352//=====================================================================
353//* CalculateExtent ---------------------------------------------------
354
355G4bool
357 const G4VoxelLimits& pVoxelLimit,
358 const G4AffineTransform& pTransform,
359 G4double& pMin, G4double& pMax) const
360{
361 G4ThreeVector bmin, bmax;
362
363 // Get bounding box
364 BoundingLimits(bmin,bmax);
365
366 // Find extent
367 G4BoundingEnvelope bbox(bmin,bmax);
368 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
369}
370
371
372//=====================================================================
373//* Inside ------------------------------------------------------------
374
376{
377
378 const G4double halftol
380 // static G4int timerid = -1;
381 // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
382 // timer.Start();
383
384 G4ThreeVector *tmpp;
385 EInside *tmpinside;
386 if (fLastInside.p == p)
387 {
388 return fLastInside.inside;
389 }
390 else
391 {
392 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
393 tmpinside = const_cast<EInside*>(&(fLastInside.inside));
394 tmpp->set(p.x(), p.y(), p.z());
395 }
396
397 EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
398 G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
399 G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
400
401 if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
402 {
403 *tmpinside = kOutside;
404 }
405 else if (outerhypearea == kSurface)
406 {
407 *tmpinside = kSurface;
408 }
409 else
410 {
411 if (distanceToOut <= halftol)
412 {
413 *tmpinside = kSurface;
414 }
415 else
416 {
417 *tmpinside = kInside;
418 }
419 }
420
421 return fLastInside.inside;
422}
423
424//=====================================================================
425//* SurfaceNormal -----------------------------------------------------
426
428{
429 //
430 // return the normal unit vector to the Hyperbolical Surface at a point
431 // p on (or nearly on) the surface
432 //
433 // Which of the three or four surfaces are we closest to?
434 //
435
436 if (fLastNormal.p == p)
437 {
438 return fLastNormal.vec;
439 }
440 auto tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
441 auto tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
442 auto tmpsurface = const_cast<G4VTwistSurface**>(fLastNormal.surface);
443 tmpp->set(p.x(), p.y(), p.z());
444
445 G4double distance = kInfinity;
446
447 G4VTwistSurface *surfaces[6];
448 surfaces[0] = fLatterTwisted;
449 surfaces[1] = fFormerTwisted;
450 surfaces[2] = fInnerHype;
451 surfaces[3] = fOuterHype;
452 surfaces[4] = fLowerEndcap;
453 surfaces[5] = fUpperEndcap;
454
455 G4ThreeVector xx;
456 G4ThreeVector bestxx;
457 G4int besti = -1;
458 for (auto i=0; i<6; ++i)
459 {
460 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
461 if (tmpdistance < distance)
462 {
463 distance = tmpdistance;
464 bestxx = xx;
465 besti = i;
466 }
467 }
468
469 tmpsurface[0] = surfaces[besti];
470 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
471
472 return fLastNormal.vec;
473}
474
475//=====================================================================
476//* DistanceToIn (p, v) -----------------------------------------------
477
479 const G4ThreeVector& v ) const
480{
481
482 // DistanceToIn (p, v):
483 // Calculate distance to surface of shape from `outside'
484 // along with the v, allowing for tolerance.
485 // The function returns kInfinity if no intersection or
486 // just grazing within tolerance.
487
488 //
489 // checking last value
490 //
491
492 G4ThreeVector* tmpp;
493 G4ThreeVector* tmpv;
494 G4double* tmpdist;
495 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
496 {
497 return fLastDistanceToIn.value;
498 }
499 else
500 {
501 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
502 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
503 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
504 tmpp->set(p.x(), p.y(), p.z());
505 tmpv->set(v.x(), v.y(), v.z());
506 }
507
508 //
509 // Calculate DistanceToIn(p,v)
510 //
511
512 EInside currentside = Inside(p);
513
514 if (currentside == kInside)
515 {
516 }
517 else
518 {
519 if (currentside == kSurface)
520 {
521 // particle is just on a boundary.
522 // If the particle is entering to the volume, return 0.
523 //
524 G4ThreeVector normal = SurfaceNormal(p);
525 if (normal*v < 0)
526 {
527 *tmpdist = 0.;
528 return fLastDistanceToInWithV.value;
529 }
530 }
531 }
532
533 // now, we can take smallest positive distance.
534
535 // Initialize
536 //
537 G4double distance = kInfinity;
538
539 // find intersections and choose nearest one.
540 //
541 G4VTwistSurface* surfaces[6];
542 surfaces[0] = fLowerEndcap;
543 surfaces[1] = fUpperEndcap;
544 surfaces[2] = fLatterTwisted;
545 surfaces[3] = fFormerTwisted;
546 surfaces[4] = fInnerHype;
547 surfaces[5] = fOuterHype;
548
549 G4ThreeVector xx;
550 G4ThreeVector bestxx;
551 for (const auto & surface : surfaces)
552 {
553 G4double tmpdistance = surface->DistanceToIn(p, v, xx);
554 if (tmpdistance < distance)
555 {
556 distance = tmpdistance;
557 bestxx = xx;
558 }
559 }
560 *tmpdist = distance;
561
562 return fLastDistanceToInWithV.value;
563}
564
565//=====================================================================
566//* DistanceToIn (p) --------------------------------------------------
567
569{
570 // DistanceToIn(p):
571 // Calculate distance to surface of shape from `outside',
572 // allowing for tolerance
573
574 //
575 // checking last value
576 //
577
578 G4ThreeVector* tmpp;
579 G4double* tmpdist;
580 if (fLastDistanceToIn.p == p)
581 {
582 return fLastDistanceToIn.value;
583 }
584 else
585 {
586 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
587 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
588 tmpp->set(p.x(), p.y(), p.z());
589 }
590
591 //
592 // Calculate DistanceToIn(p)
593 //
594
595 EInside currentside = Inside(p);
596
597 switch (currentside)
598 {
599 case (kInside) :
600 {}
601 case (kSurface) :
602 {
603 *tmpdist = 0.;
604 return fLastDistanceToIn.value;
605 }
606 case (kOutside) :
607 {
608 // Initialize
609 G4double distance = kInfinity;
610
611 // find intersections and choose nearest one.
612 G4VTwistSurface *surfaces[6];
613 surfaces[0] = fLowerEndcap;
614 surfaces[1] = fUpperEndcap;
615 surfaces[2] = fLatterTwisted;
616 surfaces[3] = fFormerTwisted;
617 surfaces[4] = fInnerHype;
618 surfaces[5] = fOuterHype;
619
620 G4ThreeVector xx;
621 G4ThreeVector bestxx;
622 for (const auto & surface : surfaces)
623 {
624 G4double tmpdistance = surface->DistanceTo(p, xx);
625 if (tmpdistance < distance)
626 {
627 distance = tmpdistance;
628 bestxx = xx;
629 }
630 }
631 *tmpdist = distance;
632 return fLastDistanceToIn.value;
633 }
634 default :
635 {
636 G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
637 FatalException, "Unknown point location!");
638 }
639 } // switch end
640
641 return kInfinity;
642}
643
644//=====================================================================
645//* DistanceToOut (p, v) ----------------------------------------------
646
648 const G4ThreeVector& v,
649 const G4bool calcNorm,
650 G4bool *validNorm,
651 G4ThreeVector *norm ) const
652{
653 // DistanceToOut (p, v):
654 // Calculate distance to surface of shape from `inside'
655 // along with the v, allowing for tolerance.
656 // The function returns kInfinity if no intersection or
657 // just grazing within tolerance.
658
659 //
660 // checking last value
661 //
662
663 G4ThreeVector* tmpp;
664 G4ThreeVector* tmpv;
665 G4double* tmpdist;
666 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
667 {
668 return fLastDistanceToOutWithV.value;
669 }
670 else
671 {
672 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
673 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
674 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
675 tmpp->set(p.x(), p.y(), p.z());
676 tmpv->set(v.x(), v.y(), v.z());
677 }
678
679 //
680 // Calculate DistanceToOut(p,v)
681 //
682
683 EInside currentside = Inside(p);
684
685 if (currentside == kOutside)
686 {
687 }
688 else
689 {
690 if (currentside == kSurface)
691 {
692 // particle is just on a boundary.
693 // If the particle is exiting from the volume, return 0.
694 //
695 G4ThreeVector normal = SurfaceNormal(p);
696 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
697 if (normal*v > 0)
698 {
699 if (calcNorm)
700 {
701 *norm = (blockedsurface->GetNormal(p, true));
702 *validNorm = blockedsurface->IsValidNorm();
703 }
704 *tmpdist = 0.;
705 return fLastDistanceToOutWithV.value;
706 }
707 }
708 }
709
710 // now, we can take smallest positive distance.
711
712 // Initialize
713 //
714 G4double distance = kInfinity;
715
716 // find intersections and choose nearest one.
717 //
718 G4VTwistSurface* surfaces[6];
719 surfaces[0] = fLatterTwisted;
720 surfaces[1] = fFormerTwisted;
721 surfaces[2] = fInnerHype;
722 surfaces[3] = fOuterHype;
723 surfaces[4] = fLowerEndcap;
724 surfaces[5] = fUpperEndcap;
725
726 G4int besti = -1;
727 G4ThreeVector xx;
728 G4ThreeVector bestxx;
729 for (auto i=0; i<6; ++i)
730 {
731 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
732 if (tmpdistance < distance)
733 {
734 distance = tmpdistance;
735 bestxx = xx;
736 besti = i;
737 }
738 }
739
740 if (calcNorm)
741 {
742 if (besti != -1)
743 {
744 *norm = (surfaces[besti]->GetNormal(p, true));
745 *validNorm = surfaces[besti]->IsValidNorm();
746 }
747 }
748
749 *tmpdist = distance;
750
751 return fLastDistanceToOutWithV.value;
752}
753
754
755//=====================================================================
756//* DistanceToOut (p) ----------------------------------------------
757
759{
760 // DistanceToOut(p):
761 // Calculate distance to surface of shape from `inside',
762 // allowing for tolerance
763
764 //
765 // checking last value
766 //
767
768 G4ThreeVector* tmpp;
769 G4double* tmpdist;
770 if (fLastDistanceToOut.p == p)
771 {
772 return fLastDistanceToOut.value;
773 }
774 else
775 {
776 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
777 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
778 tmpp->set(p.x(), p.y(), p.z());
779 }
780
781 //
782 // Calculate DistanceToOut(p)
783 //
784
785 EInside currentside = Inside(p);
786
787 switch (currentside)
788 {
789 case (kOutside) :
790 {
791 }
792 case (kSurface) :
793 {
794 *tmpdist = 0.;
795 return fLastDistanceToOut.value;
796 }
797 case (kInside) :
798 {
799 // Initialize
800 G4double distance = kInfinity;
801
802 // find intersections and choose nearest one.
803 G4VTwistSurface* surfaces[6];
804 surfaces[0] = fLatterTwisted;
805 surfaces[1] = fFormerTwisted;
806 surfaces[2] = fInnerHype;
807 surfaces[3] = fOuterHype;
808 surfaces[4] = fLowerEndcap;
809 surfaces[5] = fUpperEndcap;
810
811 G4ThreeVector xx;
812 G4ThreeVector bestxx;
813 for (const auto & surface : surfaces)
814 {
815 G4double tmpdistance = surface->DistanceTo(p, xx);
816 if (tmpdistance < distance)
817 {
818 distance = tmpdistance;
819 bestxx = xx;
820 }
821 }
822 *tmpdist = distance;
823
824 return fLastDistanceToOut.value;
825 }
826 default :
827 {
828 G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
829 FatalException, "Unknown point location!");
830 }
831 } // switch end
832
833 return 0.;
834}
835
836//=====================================================================
837//* StreamInfo --------------------------------------------------------
838
839std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
840{
841 //
842 // Stream object contents to an output stream
843 //
844 G4long oldprc = os.precision(16);
845 os << "-----------------------------------------------------------\n"
846 << " *** Dump for solid - " << GetName() << " ***\n"
847 << " ===================================================\n"
848 << " Solid type: G4TwistedTubs\n"
849 << " Parameters: \n"
850 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
851 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
852 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
853 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
854 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
855 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
856 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
857 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
858 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
859 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
860 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
861 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
862 << "-----------------------------------------------------------\n";
863 os.precision(oldprc);
864
865 return os;
866}
867
868
869//=====================================================================
870//* DiscribeYourselfTo ------------------------------------------------
871
873{
874 scene.AddSolid (*this);
875}
876
877//=====================================================================
878//* GetExtent ---------------------------------------------------------
879
881{
882 // Define the sides of the box into which the G4Tubs instance would fit.
883 //
884 G4ThreeVector pmin,pmax;
885 BoundingLimits(pmin,pmax);
886 return { pmin.x(),pmax.x(),
887 pmin.y(),pmax.y(),
888 pmin.z(),pmax.z() };
889}
890
891//=====================================================================
892//* CreatePolyhedron --------------------------------------------------
893
895{
896 // number of meshes
897 //
898 G4double absPhiTwist = std::abs(fPhiTwist);
899 G4double dA = std::max(fDPhi,absPhiTwist);
900 const G4int k =
902 const G4int n =
903 G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
904
905 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
906 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
907
908 auto ph = new G4Polyhedron;
909 typedef G4double G4double3[3];
910 typedef G4int G4int4[4];
911 auto xyz = new G4double3[nnodes]; // number of nodes
912 auto faces = new G4int4[nfaces] ; // number of faces
913 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
914 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
915 fInnerHype->GetFacets(k,n,xyz,faces,2) ;
916 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
917 fOuterHype->GetFacets(k,n,xyz,faces,4) ;
918 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
919
920 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
921
922 delete[] xyz;
923 delete[] faces;
924
925 return ph;
926}
927
928//=====================================================================
929//* GetPolyhedron -----------------------------------------------------
930
932{
933 if (fpPolyhedron == nullptr ||
934 fRebuildPolyhedron ||
936 fpPolyhedron->GetNumberOfRotationSteps())
937 {
938 G4AutoLock l(&polyhedronMutex);
939 delete fpPolyhedron;
940 fpPolyhedron = CreatePolyhedron();
941 fRebuildPolyhedron = false;
942 l.unlock();
943 }
944 return fpPolyhedron;
945}
946
947//=====================================================================
948//* CreateSurfaces ----------------------------------------------------
949
950void G4TwistedTubs::CreateSurfaces()
951{
952 // create 6 surfaces of TwistedTub
953
954 G4ThreeVector x0(0, 0, fEndZ[0]);
955 G4ThreeVector n (0, 0, -1);
956
957 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
958 fEndInnerRadius, fEndOuterRadius,
959 fDPhi, fEndPhi, fEndZ, -1) ;
960
961 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
962 fEndInnerRadius, fEndOuterRadius,
963 fDPhi, fEndPhi, fEndZ, 1) ;
964
965 G4RotationMatrix rotHalfDPhi;
966 rotHalfDPhi.rotateZ(0.5*fDPhi);
967
968 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
969 fEndInnerRadius, fEndOuterRadius,
970 fDPhi, fEndPhi, fEndZ,
971 fInnerRadius, fOuterRadius, fKappa,
972 1 ) ;
973 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
974 fEndInnerRadius, fEndOuterRadius,
975 fDPhi, fEndPhi, fEndZ,
976 fInnerRadius, fOuterRadius, fKappa,
977 -1 ) ;
978
979 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
980 fEndInnerRadius, fEndOuterRadius,
981 fDPhi, fEndPhi, fEndZ,
982 fInnerRadius, fOuterRadius,fKappa,
983 fTanInnerStereo, fTanOuterStereo, -1) ;
984 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
985 fEndInnerRadius, fEndOuterRadius,
986 fDPhi, fEndPhi, fEndZ,
987 fInnerRadius, fOuterRadius,fKappa,
988 fTanInnerStereo, fTanOuterStereo, 1) ;
989
990
991 // set neighbour surfaces
992 //
993 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
994 fOuterHype, fFormerTwisted);
995 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
996 fOuterHype, fFormerTwisted);
997 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
998 fOuterHype, fUpperEndcap);
999 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1000 fOuterHype, fUpperEndcap);
1001 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1002 fFormerTwisted, fUpperEndcap);
1003 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1004 fFormerTwisted, fUpperEndcap);
1005}
1006
1007
1008//=====================================================================
1009//* GetEntityType -----------------------------------------------------
1010
1012{
1013 return {"G4TwistedTubs"};
1014}
1015
1016//=====================================================================
1017//* Clone -------------------------------------------------------------
1018
1020{
1021 return new G4TwistedTubs(*this);
1022}
1023
1024//=====================================================================
1025//* GetCubicVolume ----------------------------------------------------
1026
1028{
1029 if (fCubicVolume == 0.)
1030 {
1031 G4double DPhi = GetDPhi();
1032 G4double Z0 = GetEndZ(0);
1033 G4double Z1 = GetEndZ(1);
1034 G4double Ain = GetInnerRadius();
1035 G4double Aout = GetOuterRadius();
1036 G4double R0in = GetEndInnerRadius(0);
1037 G4double R1in = GetEndInnerRadius(1);
1038 G4double R0out = GetEndOuterRadius(0);
1039 G4double R1out = GetEndOuterRadius(1);
1040
1041 // V_hyperboloid = pi*h*(2*a*a + R*R)/3
1042 fCubicVolume = (2.*(Z1 - Z0)*(Aout + Ain)*(Aout - Ain)
1043 + Z1*(R1out + R1in)*(R1out - R1in)
1044 - Z0*(R0out + R0in)*(R0out - R0in))*DPhi/6.;
1045 }
1046 return fCubicVolume;
1047}
1048
1049//=====================================================================
1050//* GetLateralArea ----------------------------------------------------
1051
1053G4TwistedTubs::GetLateralArea(G4double a, G4double r, G4double z) const
1054{
1055 if (z == 0) return 0.;
1056 G4double h = std::abs(z);
1057 G4double area = h*a;
1058 if (std::abs(a - r) > kCarTolerance)
1059 {
1060 G4double aa = a*a;
1061 G4double hh = h*h;
1062 G4double rr = r*r;
1063 G4double cc = aa*hh/(rr - aa);
1064 G4double k = std::sqrt(aa + cc)/cc;
1065 G4double kh = k*h;
1066 area = 0.5*a*(h*std::sqrt(1. + kh*kh) + std::asinh(kh)/k);
1067 }
1068 return GetDPhi()*area;
1069}
1070
1071//=====================================================================
1072//* GetPhiCutArea -----------------------------------------------------
1073
1075G4TwistedTubs::GetPhiCutArea(G4double a, G4double r, G4double z) const
1076{
1077 if (GetDPhi() >= CLHEP::twopi || r <= 0 || z == 0) return 0.;
1078 G4double h = std::abs(z);
1079 G4double area = h*a;
1080 if (GetPhiTwist() > kCarTolerance)
1081 {
1082 G4double sinw = std::sin(0.5*GetPhiTwist())*h/GetZHalfLength();
1083 G4double p = sinw*r/h;
1084 G4double q = sinw*r/a;
1085 G4double pp = p*p;
1086 G4double qq = q*q;
1087 G4double pq = p*q;
1088 G4double sqroot = std::sqrt(pp + qq + 1);
1089 area = (pq*sqroot +
1090 0.5*p*(pp + 3.)*std::atanh(q/sqroot) +
1091 0.5*q*(qq + 3.)*std::atanh(p/sqroot) +
1092 std::atan(sqroot/(pq)) - CLHEP::halfpi)*h*a/(3.*pq);
1093 }
1094 return area;
1095}
1096
1097//=====================================================================
1098//* GetSurfaceArea ----------------------------------------------------
1099
1101{
1102 if (fSurfaceArea == 0.)
1103 {
1104 G4double dphi = GetDPhi();
1105 G4double Ainn = GetInnerRadius();
1106 G4double Aout = GetOuterRadius();
1107 G4double Rinn0 = GetEndInnerRadius(0);
1108 G4double Rout0 = GetEndOuterRadius(0);
1109 G4double Rinn1 = GetEndInnerRadius(1);
1110 G4double Rout1 = GetEndOuterRadius(1);
1111 G4double z0 = GetEndZ(0);
1112 G4double z1 = GetEndZ(1);
1113
1114 G4double base0 = 0.5*dphi*(Rout0*Rout0 - Rinn0*Rinn0); // lower base
1115 G4double inner0 = GetLateralArea(Ainn, Rinn0, z0); // lower inner surface
1116 G4double outer0 = GetLateralArea(Aout, Rout0, z0); // lower outer surface
1117 G4double cut0 = // lower phi cut
1118 GetPhiCutArea(Aout, Rout0, z0) - GetPhiCutArea(Ainn, Rinn0, z0);
1119
1120 G4double base1 = base0;
1121 G4double inner1 = inner0;
1122 G4double outer1 = outer0;
1123 G4double cut1 = cut0;
1124 if (std::abs(z0) != std::abs(z1))
1125 {
1126 base1 = 0.5*dphi*(Rout1*Rout1 - Rinn1*Rinn1); // upper base
1127 inner1 = GetLateralArea(Ainn, Rinn1, z1); // upper inner surface
1128 outer1 = GetLateralArea(Aout, Rout1, z1); // upper outer surface
1129 cut1 = // upper phi cut
1130 GetPhiCutArea(Aout, Rout1, z1) - GetPhiCutArea(Ainn, Rinn1, z1);
1131 }
1132 fSurfaceArea = base0 + base1 +
1133 ((z0*z1 < 0) ?
1134 (inner0 + inner1 + outer0 + outer1 + 2.*(cut0 + cut1)) :
1135 std::abs(inner0 - inner1 + outer0 - outer1 + 2.*(cut0 - cut1)));
1136 }
1137 return fSurfaceArea;
1138}
1139
1140//=====================================================================
1141//* GetPointOnSurface -------------------------------------------------
1142
1144{
1145
1146 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1147 G4double phi , phimin, phimax ;
1148 G4double x , xmin, xmax ;
1149 G4double r , rmin, rmax ;
1150
1151 G4double a1 = fOuterHype->GetSurfaceArea() ;
1152 G4double a2 = fInnerHype->GetSurfaceArea() ;
1153 G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1154 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1155 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1156 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1157
1158 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1159
1160 if(chose < a1)
1161 {
1162
1163 phimin = fOuterHype->GetBoundaryMin(z) ;
1164 phimax = fOuterHype->GetBoundaryMax(z) ;
1165 phi = G4RandFlat::shoot(phimin,phimax) ;
1166
1167 return fOuterHype->SurfacePoint(phi,z,true) ;
1168
1169 }
1170 else if ( (chose >= a1) && (chose < a1 + a2 ) )
1171 {
1172
1173 phimin = fInnerHype->GetBoundaryMin(z) ;
1174 phimax = fInnerHype->GetBoundaryMax(z) ;
1175 phi = G4RandFlat::shoot(phimin,phimax) ;
1176
1177 return fInnerHype->SurfacePoint(phi,z,true) ;
1178
1179 }
1180 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1181 {
1182
1183 xmin = fLatterTwisted->GetBoundaryMin(z) ;
1184 xmax = fLatterTwisted->GetBoundaryMax(z) ;
1185 x = G4RandFlat::shoot(xmin,xmax) ;
1186
1187 return fLatterTwisted->SurfacePoint(x,z,true) ;
1188
1189 }
1190 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1191 {
1192
1193 xmin = fFormerTwisted->GetBoundaryMin(z) ;
1194 xmax = fFormerTwisted->GetBoundaryMax(z) ;
1195 x = G4RandFlat::shoot(xmin,xmax) ;
1196
1197 return fFormerTwisted->SurfacePoint(x,z,true) ;
1198 }
1199 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1200 {
1201 rmin = GetEndInnerRadius(0) ;
1202 rmax = GetEndOuterRadius(0) ;
1203 r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1204
1205 phimin = fLowerEndcap->GetBoundaryMin(r) ;
1206 phimax = fLowerEndcap->GetBoundaryMax(r) ;
1207 phi = G4RandFlat::shoot(phimin,phimax) ;
1208
1209 return fLowerEndcap->SurfacePoint(phi,r,true) ;
1210 }
1211 else
1212 {
1213 rmin = GetEndInnerRadius(1) ;
1214 rmax = GetEndOuterRadius(1) ;
1215 r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1216
1217 phimin = fUpperEndcap->GetBoundaryMin(r) ;
1218 phimax = fUpperEndcap->GetBoundaryMax(r) ;
1219 phi = G4RandFlat::shoot(phimin,phimax) ;
1220
1221 return fUpperEndcap->SurfacePoint(phi,r,true) ;
1222 }
1223}
@ JustWarning
@ FatalException
@ FatalErrorInArgument
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
#define G4MUTEX_INITIALIZER
std::mutex G4Mutex
double G4double
Definition G4Types.hh:83
long G4long
Definition G4Types.hh:87
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
double x() const
double y() const
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
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
G4double GetRadialTolerance() const
static G4GeometryTolerance * GetInstance()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=nullptr, G4ThreeVector *n=nullptr) const override
G4double GetOuterRadius() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const override
G4double GetZHalfLength() const
G4double GetPhiTwist() const
G4Polyhedron * CreatePolyhedron() const override
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const override
G4TwistedTubs & operator=(const G4TwistedTubs &rhs)
G4double GetEndInnerRadius() const
G4double GetEndOuterRadius() const
G4VSolid * Clone() const override
G4Polyhedron * GetPolyhedron() const override
~G4TwistedTubs() override
G4GeometryType GetEntityType() const override
G4double GetEndPhi(G4int i) const
G4TwistedTubs(const G4String &pname, G4double twistedangle, G4double endinnerrad, G4double endouterrad, G4double halfzlen, G4double dphi)
std::ostream & StreamInfo(std::ostream &os) const override
G4double GetSurfaceArea() override
G4double GetCubicVolume() override
void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *) override
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
G4VisExtent GetExtent() const override
G4ThreeVector GetPointOnSurface() const override
void DescribeYourselfTo(G4VGraphicsScene &scene) const override
G4double GetEndZ(G4int i) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const override
G4double GetInnerRadius() const
EInside Inside(const G4ThreeVector &p) const override
G4double GetDPhi() const
virtual void AddSolid(const G4Box &)=0
G4String GetName() const
void DumpInfo() const
G4double kCarTolerance
Definition G4VSolid.hh:299
G4VSolid & operator=(const G4VSolid &rhs)
Definition G4VSolid.cc:107
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 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()
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