| 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 // G4VTwistedFaceted implementation
27 //
28 // Author: 04-Nov-2004 - O.Link (Oliver.Link@c
29 // -------------------------------------------
30
31 #include "G4VTwistedFaceted.hh"
32
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4VoxelLimits.hh"
36 #include "G4AffineTransform.hh"
37 #include "G4BoundingEnvelope.hh"
38 #include "G4SolidExtentList.hh"
39 #include "G4ClippablePolygon.hh"
40 #include "G4VPVParameterisation.hh"
41 #include "G4GeometryTolerance.hh"
42 #include "meshdefs.hh"
43
44 #include "G4VGraphicsScene.hh"
45 #include "G4Polyhedron.hh"
46 #include "G4VisExtent.hh"
47
48 #include "Randomize.hh"
49
50 #include "G4AutoLock.hh"
51
52 namespace
53 {
54 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
55 }
56
57
58 //============================================
59 //* constructors -----------------------------
60
61 G4VTwistedFaceted::
62 G4VTwistedFaceted( const G4String& pname,
63 G4double PhiTwist,
64 G4double pDz,
65 G4double pTheta, //
66 G4double pPhi, //
67 G4double pDy1, //
68 G4double pDx1, //
69 G4double pDx2, //
70 G4double pDy2, //
71 G4double pDx3, //
72 G4double pDx4, //
73 G4double pAlph ) //
74 : G4VSolid(pname),
75 fLowerEndcap(nullptr), fUpperEndcap(nullpt
76 fSide90(nullptr), fSide180(nullptr), fSide
77 {
78
79 G4double pDytmp ;
80 G4double fDxUp ;
81 G4double fDxDown ;
82
83 fDx1 = pDx1 ;
84 fDx2 = pDx2 ;
85 fDx3 = pDx3 ;
86 fDx4 = pDx4 ;
87 fDy1 = pDy1 ;
88 fDy2 = pDy2 ;
89 fDz = pDz ;
90
91 G4double kAngTolerance
92 = G4GeometryTolerance::GetInstance()->GetA
93
94 // maximum values
95 //
96 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
97 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
98 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDow
99 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
100
101 // planarity check
102 //
103 if ( fDx1 != fDx2 && fDx3 != fDx4 )
104 {
105 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 -
106 if ( std::fabs(pDytmp - fDy2) > kCarTolera
107 {
108 std::ostringstream message;
109 message << "Not planar surface in untwis
110 << GetName() << G4endl
111 << "fDy2 is " << fDy2 << " but s
112 << pDytmp << ".";
113 G4Exception("G4VTwistedFaceted::G4VTwist
114 FatalErrorInArgument, messag
115 }
116 }
117
118 #ifdef G4TWISTDEBUG
119 if ( fDx1 == fDx2 && fDx3 == fDx4 )
120 {
121 G4cout << "Trapezoid is a box" << G4endl
122 }
123
124 #endif
125
126 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || (
127 {
128 std::ostringstream message;
129 message << "Not planar surface in untwiste
130 << GetName() << G4endl
131 << "One endcap is rectangular, the
132 << "For planarity reasons they hav
133 << "on both sides.";
134 G4Exception("G4VTwistedFaceted::G4VTwisted
135 FatalErrorInArgument, message)
136 }
137
138 // twist angle
139 //
140 fPhiTwist = PhiTwist ;
141
142 // tilt angle
143 //
144 fAlph = pAlph ;
145 fTAlph = std::tan(fAlph) ;
146
147 fTheta = pTheta ;
148 fPhi = pPhi ;
149
150 // dx in surface equation
151 //
152 fdeltaX = 2 * fDz * std::tan(fTheta) * std::
153
154 // dy in surface equation
155 //
156 fdeltaY = 2 * fDz * std::tan(fTheta) * std::
157
158 if ( ( fDx1 <= 2*kCarTolerance)
159 || ( fDx2 <= 2*kCarTolerance)
160 || ( fDx3 <= 2*kCarTolerance)
161 || ( fDx4 <= 2*kCarTolerance)
162 || ( fDy1 <= 2*kCarTolerance)
163 || ( fDy2 <= 2*kCarTolerance)
164 || ( fDz <= 2*kCarTolerance)
165 || ( std::fabs(fPhiTwist) <= 2*kAngTo
166 || ( std::fabs(fPhiTwist) >= pi/2 )
167 || ( std::fabs(fAlph) >= pi/2 )
168 || fTheta >= pi/2 || fTheta < 0
169 )
170 {
171 std::ostringstream message;
172 message << "Invalid dimensions. Too small,
173 << GetName() << G4endl
174 << "fDx 1-4 = " << fDx1/cm << ", "
175 << fDx3/cm << ", " << fDx4/cm << "
176 << "fDy 1-2 = " << fDy1/cm << ", "
177 << " cm" << G4endl
178 << "fDz = " << fDz/cm << " cm" <<
179 << " twistangle " << fPhiTwist/deg
180 << " phi,theta = " << fPhi/deg <<
181 G4Exception("G4TwistedTrap::G4VTwistedFace
182 "GeomSolids0002", FatalErrorIn
183 }
184 CreateSurfaces();
185 }
186
187
188 //============================================
189 //* Fake default constructor -----------------
190
191 G4VTwistedFaceted::G4VTwistedFaceted( __void__
192 : G4VSolid(a),
193 fLowerEndcap(nullptr), fUpperEndcap(nullpt
194 fSide0(nullptr), fSide90(nullptr), fSide18
195 {
196 }
197
198
199 //============================================
200 //* destructor -------------------------------
201
202 G4VTwistedFaceted::~G4VTwistedFaceted()
203 {
204 delete fLowerEndcap ;
205 delete fUpperEndcap ;
206
207 delete fSide0 ;
208 delete fSide90 ;
209 delete fSide180 ;
210 delete fSide270 ;
211 delete fpPolyhedron; fpPolyhedron = nullptr;
212 }
213
214
215 //============================================
216 //* Copy constructor -------------------------
217
218 G4VTwistedFaceted::G4VTwistedFaceted(const G4V
219 : G4VSolid(rhs),
220 fCubicVolume(rhs.fCubicVolume), fSurfaceAr
221 fTheta(rhs.fTheta), fPhi(rhs.fPhi),
222 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f
223 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fD
224 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdel
225 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTw
226 fUpperEndcap(nullptr), fSide0(nullptr), fS
227 {
228 CreateSurfaces();
229 }
230
231
232 //============================================
233 //* Assignment operator ----------------------
234
235 G4VTwistedFaceted& G4VTwistedFaceted::operator
236 {
237 // Check assignment to self
238 //
239 if (this == &rhs) { return *this; }
240
241 // Copy base class data
242 //
243 G4VSolid::operator=(rhs);
244
245 // Copy data
246 //
247 fTheta = rhs.fTheta; fPhi = rhs.fPhi;
248 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.f
249 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fD
250 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdelt
251 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTw
252 fUpperEndcap= nullptr; fSide0= nullptr; fSi
253 fCubicVolume= rhs.fCubicVolume; fSurfaceAre
254 fRebuildPolyhedron = false;
255 delete fpPolyhedron; fpPolyhedron = nullptr
256
257 CreateSurfaces();
258
259 return *this;
260 }
261
262
263 //============================================
264 //* ComputeDimensions ------------------------
265
266 void G4VTwistedFaceted::ComputeDimensions(G4VP
267 cons
268 cons
269 {
270 G4Exception("G4VTwistedFaceted::ComputeDimen
271 "GeomSolids0001", FatalException
272 "G4VTwistedFaceted does not supp
273 }
274
275
276 //============================================
277 //* Extent -----------------------------------
278
279 void G4VTwistedFaceted::BoundingLimits(G4Three
280 G4Three
281 {
282 G4double cosPhi = std::cos(fPhi);
283 G4double sinPhi = std::sin(fPhi);
284 G4double tanTheta = std::tan(fTheta);
285 G4double tanAlpha = fTAlph;
286
287 G4double xmid1 = fDy1*tanAlpha;
288 G4double x1 = std::abs(xmid1 + fDx1);
289 G4double x2 = std::abs(xmid1 - fDx1);
290 G4double x3 = std::abs(xmid1 + fDx2);
291 G4double x4 = std::abs(xmid1 - fDx2);
292 G4double xmax1 = std::max(std::max(std::max(
293 G4double rmax1 = std::sqrt(xmax1*xmax1 + fDy
294
295 G4double xmid2 = fDy2*tanAlpha;
296 G4double x5 = std::abs(xmid2 + fDx3);
297 G4double x6 = std::abs(xmid2 - fDx3);
298 G4double x7 = std::abs(xmid2 + fDx4);
299 G4double x8 = std::abs(xmid2 - fDx4);
300 G4double xmax2 = std::max(std::max(std::max(
301 G4double rmax2 = std::sqrt(xmax2*xmax2 + fDy
302
303 G4double x0 = fDz*tanTheta*cosPhi;
304 G4double y0 = fDz*tanTheta*sinPhi;
305 G4double xmin = std::min(-x0 - rmax1, x0 - r
306 G4double ymin = std::min(-y0 - rmax1, y0 - r
307 G4double xmax = std::max(-x0 + rmax1, x0 + r
308 G4double ymax = std::max(-y0 + rmax1, y0 + r
309 pMin.set(xmin, ymin,-fDz);
310 pMax.set(xmax, ymax, fDz);
311 }
312
313
314 //============================================
315 //* CalculateExtent --------------------------
316
317 G4bool
318 G4VTwistedFaceted::CalculateExtent( const EAxi
319 const G4Vo
320 const G4Af
321 G4do
322 G4do
323 {
324 G4ThreeVector bmin, bmax;
325
326 // Get bounding box
327 BoundingLimits(bmin,bmax);
328
329 // Find extent
330 G4BoundingEnvelope bbox(bmin,bmax);
331 return bbox.CalculateExtent(pAxis,pVoxelLimi
332 }
333
334
335 //============================================
336 //* Inside -----------------------------------
337
338 EInside G4VTwistedFaceted::Inside(const G4Thre
339 {
340
341 EInside tmpin = kOutside ;
342
343 G4double phi = p.z()/(2*fDz) * fPhiTwist ;
344 G4double cphi = std::cos(-phi) ;
345 G4double sphi = std::sin(-phi) ;
346
347 G4double px = p.x() + fdeltaX * ( -phi/fPh
348 G4double py = p.y() + fdeltaY * ( -phi/fPh
349 G4double pz = p.z() ;
350
351 G4double posx = px * cphi - py * sphi ;
352 G4double posy = px * sphi + py * cphi ;
353 G4double posz = pz ;
354
355 G4double xMin = Xcoef(posy,phi,fTAlph) - 2*
356 G4double xMax = Xcoef(posy,phi,fTAlph) ;
357
358 G4double yMax = GetValueB(phi)/2. ; // b(p
359 G4double yMin = -yMax ;
360
361 #ifdef G4TWISTDEBUG
362
363 G4cout << "inside called: p = " << p << G4e
364 G4cout << "fDx1 = " << fDx1 << G4endl ;
365 G4cout << "fDx2 = " << fDx2 << G4endl ;
366 G4cout << "fDx3 = " << fDx3 << G4endl ;
367 G4cout << "fDx4 = " << fDx4 << G4endl ;
368
369 G4cout << "fDy1 = " << fDy1 << G4endl ;
370 G4cout << "fDy2 = " << fDy2 << G4endl ;
371
372 G4cout << "fDz = " << fDz << G4endl ;
373
374 G4cout << "Tilt angle alpha = " << fAlph <<
375 G4cout << "phi,theta = " << fPhi << " , " <
376
377 G4cout << "Twist angle = " << fPhiTwist <<
378
379 G4cout << "posx = " << posx << G4endl ;
380 G4cout << "posy = " << posy << G4endl ;
381 G4cout << "xMin = " << xMin << G4endl ;
382 G4cout << "xMax = " << xMax << G4endl ;
383 G4cout << "yMin = " << yMin << G4endl ;
384 G4cout << "yMax = " << yMax << G4endl ;
385
386 #endif
387
388
389 if ( posx <= xMax - kCarTolerance*0.5
390 && posx >= xMin + kCarTolerance*0.5 )
391 {
392 if ( posy <= yMax - kCarTolerance*0.5
393 && posy >= yMin + kCarTolerance*0.5 )
394 {
395 if (std::fabs(posz) <= fDz - kCarTo
396 else if (std::fabs(posz) <= fDz + kCarTo
397 }
398 else if ( posy <= yMax + kCarTolerance*0.5
399 && posy >= yMin - kCarTolerance*0.5
400 {
401 if (std::fabs(posz) <= fDz + kCarToleran
402 }
403 }
404 else if ( posx <= xMax + kCarTolerance*0.5
405 && posx >= xMin - kCarTolerance*0.5 )
406 {
407 if ( posy <= yMax + kCarTolerance*0.5
408 && posy >= yMin - kCarTolerance*0.5 )
409 {
410 if (std::fabs(posz) <= fDz + kCarToleran
411 }
412 }
413
414 #ifdef G4TWISTDEBUG
415 G4cout << "inside = " << tmpin << G4endl ;
416 #endif
417
418 return tmpin;
419
420 }
421
422
423 //============================================
424 //* SurfaceNormal ----------------------------
425
426 G4ThreeVector G4VTwistedFaceted::SurfaceNormal
427 {
428 //
429 // return the normal unit vector to the Hyp
430 // p on (or nearly on) the surface
431 //
432 // Which of the three or four surfaces are
433 //
434
435
436 G4double distance = kInfinity;
437
438 G4VTwistSurface* surfaces[6];
439
440 surfaces[0] = fSide0 ;
441 surfaces[1] = fSide90 ;
442 surfaces[2] = fSide180 ;
443 surfaces[3] = fSide270 ;
444 surfaces[4] = fLowerEndcap;
445 surfaces[5] = fUpperEndcap;
446
447 G4ThreeVector xx;
448 G4ThreeVector bestxx;
449 G4int i;
450 G4int besti = -1;
451 for (i=0; i< 6; i++)
452 {
453 G4double tmpdistance = surfaces[i]->Dist
454 if (tmpdistance < distance)
455 {
456 distance = tmpdistance;
457 bestxx = xx;
458 besti = i;
459 }
460 }
461
462 return surfaces[besti]->GetNormal(bestxx, t
463 }
464
465
466 //============================================
467 //* DistanceToIn (p, v) ----------------------
468
469 G4double G4VTwistedFaceted::DistanceToIn (cons
470 cons
471 {
472
473 // DistanceToIn (p, v):
474 // Calculate distance to surface of shape f
475 // along with the v, allowing for tolerance
476 // The function returns kInfinity if no int
477 // just grazing within tolerance.
478
479 //
480 // Calculate DistanceToIn(p,v)
481 //
482
483 EInside currentside = Inside(p);
484
485 if (currentside == kInside)
486 {
487 }
488 else if (currentside == kSurface)
489 {
490 // particle is just on a boundary.
491 // if the particle is entering to the vol
492 //
493 G4ThreeVector normal = SurfaceNormal(p);
494 if (normal*v < 0)
495 {
496 return 0;
497 }
498 }
499
500 // now, we can take smallest positive dista
501
502 // Initialize
503 //
504 G4double distance = kInfinity;
505
506 // Find intersections and choose nearest on
507 //
508 G4VTwistSurface *surfaces[6];
509
510 surfaces[0] = fSide0;
511 surfaces[1] = fSide90 ;
512 surfaces[2] = fSide180 ;
513 surfaces[3] = fSide270 ;
514 surfaces[4] = fLowerEndcap;
515 surfaces[5] = fUpperEndcap;
516
517 G4ThreeVector xx;
518 G4ThreeVector bestxx;
519 for (const auto & surface : surfaces)
520 {
521 #ifdef G4TWISTDEBUG
522 G4cout << G4endl << "surface " << &surfa
523 #endif
524 G4double tmpdistance = surface->Distance
525 #ifdef G4TWISTDEBUG
526 G4cout << "Solid DistanceToIn : distance
527 G4cout << "intersection point = " << xx
528 #endif
529 if (tmpdistance < distance)
530 {
531 distance = tmpdistance;
532 bestxx = xx;
533 }
534 }
535
536 #ifdef G4TWISTDEBUG
537 G4cout << "best distance = " << distance <<
538 #endif
539
540 // timer.Stop();
541 return distance;
542 }
543
544
545 //============================================
546 //* DistanceToIn (p) -------------------------
547
548 G4double G4VTwistedFaceted::DistanceToIn (cons
549 {
550 // DistanceToIn(p):
551 // Calculate distance to surface of shape f
552 // allowing for tolerance
553 //
554
555 //
556 // Calculate DistanceToIn(p)
557 //
558
559 EInside currentside = Inside(p);
560
561 switch (currentside)
562 {
563 case (kInside) :
564 {
565 }
566
567 case (kSurface) :
568 {
569 return 0;
570 }
571
572 case (kOutside) :
573 {
574 // Initialize
575 //
576 G4double distance = kInfinity;
577
578 // Find intersections and choose near
579 //
580 G4VTwistSurface* surfaces[6];
581
582 surfaces[0] = fSide0;
583 surfaces[1] = fSide90 ;
584 surfaces[2] = fSide180 ;
585 surfaces[3] = fSide270 ;
586 surfaces[4] = fLowerEndcap;
587 surfaces[5] = fUpperEndcap;
588
589 G4ThreeVector xx;
590 G4ThreeVector bestxx;
591 for (const auto & surface : surfaces)
592 {
593 G4double tmpdistance = surface->Di
594 if (tmpdistance < distance)
595 {
596 distance = tmpdistance;
597 bestxx = xx;
598 }
599 }
600 return distance;
601 }
602
603 default:
604 {
605 G4Exception("G4VTwistedFaceted::Dista
606 FatalException, "Unknown
607 }
608 } // switch end
609
610 return 0.;
611 }
612
613
614 //============================================
615 //* DistanceToOut (p, v) ---------------------
616
617 G4double
618 G4VTwistedFaceted::DistanceToOut( const G4Thre
619 const G4Thre
620 const G4bool
621 G4bool
622 G4Thre
623 {
624 // DistanceToOut (p, v):
625 // Calculate distance to surface of shape f
626 // along with the v, allowing for tolerance
627 // The function returns kInfinity if no int
628 // just grazing within tolerance.
629
630 //
631 // Calculate DistanceToOut(p,v)
632 //
633
634 EInside currentside = Inside(p);
635
636 if (currentside == kOutside)
637 {
638 }
639 else if (currentside == kSurface)
640 {
641 // particle is just on a boundary.
642 // if the particle is exiting from the v
643 //
644 G4ThreeVector normal = SurfaceNormal(p);
645 if (normal*v > 0)
646 {
647 if (calcNorm)
648 {
649 *norm = normal;
650 *validNorm = true;
651 }
652 // timer.Stop();
653 return 0;
654 }
655 }
656
657 // now, we can take smallest positive dista
658
659 // Initialize
660 G4double distance = kInfinity;
661
662 // find intersections and choose nearest on
663 G4VTwistSurface *surfaces[6];
664
665 surfaces[0] = fSide0;
666 surfaces[1] = fSide90 ;
667 surfaces[2] = fSide180 ;
668 surfaces[3] = fSide270 ;
669 surfaces[4] = fLowerEndcap;
670 surfaces[5] = fUpperEndcap;
671
672 G4int besti = -1;
673 G4ThreeVector xx;
674 G4ThreeVector bestxx;
675 for (auto i=0; i<6 ; ++i)
676 {
677 G4double tmpdistance = surfaces[i]->Dist
678 if (tmpdistance < distance)
679 {
680 distance = tmpdistance;
681 bestxx = xx;
682 besti = i;
683 }
684 }
685
686 if (calcNorm)
687 {
688 if (besti != -1)
689 {
690 *norm = (surfaces[besti]->GetNormal(p
691 *validNorm = surfaces[besti]->IsValid
692 }
693 }
694
695 return distance;
696 }
697
698
699 //============================================
700 //* DistanceToOut (p) ------------------------
701
702 G4double G4VTwistedFaceted::DistanceToOut( con
703 {
704 // DistanceToOut(p):
705 // Calculate distance to surface of shape f
706 // allowing for tolerance
707
708 //
709 // Calculate DistanceToOut(p)
710 //
711
712 EInside currentside = Inside(p);
713 G4double retval = kInfinity;
714
715 switch (currentside)
716 {
717 case (kOutside) :
718 {
719 #ifdef G4SPECSDEBUG
720 G4int oldprc = G4cout.precision(16) ;
721 G4cout << G4endl ;
722 DumpInfo();
723 G4cout << "Position:" << G4endl << G4
724 G4cout << "p.x() = " << p.x()/mm <<
725 G4cout << "p.y() = " << p.y()/mm <<
726 G4cout << "p.z() = " << p.z()/mm <<
727 G4cout.precision(oldprc) ;
728 G4Exception("G4VTwistedFaceted::Distan
729 JustWarning, "Point p is o
730 #endif
731 break;
732 }
733 case (kSurface) :
734 {
735 retval = 0;
736 break;
737 }
738
739 case (kInside) :
740 {
741 // Initialize
742 //
743 G4double distance = kInfinity;
744
745 // find intersections and choose neare
746 //
747 G4VTwistSurface* surfaces[6];
748
749 surfaces[0] = fSide0;
750 surfaces[1] = fSide90 ;
751 surfaces[2] = fSide180 ;
752 surfaces[3] = fSide270 ;
753 surfaces[4] = fLowerEndcap;
754 surfaces[5] = fUpperEndcap;
755
756 G4ThreeVector xx;
757 G4ThreeVector bestxx;
758 for (const auto & surface : surfaces)
759 {
760 G4double tmpdistance = surface->Dist
761 if (tmpdistance < distance)
762 {
763 distance = tmpdistance;
764 bestxx = xx;
765 }
766 }
767 retval = distance;
768 break;
769 }
770
771 default :
772 {
773 G4Exception("G4VTwistedFaceted::Distan
774 FatalException, "Unknown p
775 break;
776 }
777 } // switch end
778
779 return retval;
780 }
781
782
783 //============================================
784 //* StreamInfo -------------------------------
785
786 std::ostream& G4VTwistedFaceted::StreamInfo(st
787 {
788 //
789 // Stream object contents to an output strea
790 //
791 G4long oldprc = os.precision(16);
792 os << "-------------------------------------
793 << " *** Dump for solid - " << GetName
794 << " =================================
795 << " Solid type: G4VTwistedFaceted\n"
796 << " Parameters: \n"
797 << " polar angle theta = " << fTheta/
798 << " azimuthal angle phi = " << fPhi/de
799 << " tilt angle alpha = " << fAlph/de
800 << " TWIST angle = " << fPhiTwis
801 << " Half length along y (lower endcap)
802 << G4endl
803 << " Half length along x (lower endcap,
804 << G4endl
805 << " Half length along x (lower endcap,
806 << G4endl
807 << " Half length along y (upper endcap)
808 << G4endl
809 << " Half length along x (upper endcap,
810 << G4endl
811 << " Half length along x (upper endcap,
812 << G4endl
813 << "-------------------------------------
814 os.precision(oldprc);
815
816 return os;
817 }
818
819
820 //============================================
821 //* DiscribeYourselfTo -----------------------
822
823 void G4VTwistedFaceted::DescribeYourselfTo (G4
824 {
825 scene.AddSolid (*this);
826 }
827
828
829 //============================================
830 //* GetExtent --------------------------------
831
832 G4VisExtent G4VTwistedFaceted::GetExtent() con
833 {
834 G4double maxRad = std::sqrt( fDx*fDx + fDy*f
835
836 return { -maxRad, maxRad ,
837 -maxRad, maxRad ,
838 -fDz, fDz };
839 }
840
841
842 //============================================
843 //* CreateSurfaces ---------------------------
844
845 void G4VTwistedFaceted::CreateSurfaces()
846 {
847
848 // create 6 surfaces of TwistedTub.
849
850 if ( fDx1 == fDx2 && fDx3 == fDx4 ) // sp
851 {
852 fSide0 = new G4TwistBoxSide("0deg", fPhi
853 fDy1, fDx1, fDx1, fD
854 fSide180 = new G4TwistBoxSide("180deg", fP
855 fDy1, fDx1, fDx1, fD
856 }
857 else // default general case
858 {
859 fSide0 = new G4TwistTrapAlphaSide("0deg"
860 fPhi, fDy1, fDx1, fDx2,
861 fSide180 = new G4TwistTrapAlphaSide("180de
862 fPhi+pi, fDy1, fDx2, fDx1, fD
863 }
864
865 // create parallel sides
866 //
867 fSide90 = new G4TwistTrapParallelSide("90deg
868 fPhi, fDy1, fDx1, fDx2,
869 fSide270 = new G4TwistTrapParallelSide("270d
870 fPhi+pi, fDy1, fDx2, fDx1, fD
871
872 // create endcaps
873 //
874 fUpperEndcap = new G4TwistTrapFlatSide("Uppe
875 fDz, fAlph
876 fLowerEndcap = new G4TwistTrapFlatSide("Lowe
877 fDz, fAlph
878
879 // Set neighbour surfaces
880
881 fSide0->SetNeighbours( fSide270 , fLowerEnd
882 fSide90->SetNeighbours( fSide0 , fLowerEnd
883 fSide180->SetNeighbours(fSide90 , fLowerEnd
884 fSide270->SetNeighbours(fSide180 , fLowerEnd
885 fUpperEndcap->SetNeighbours( fSide180, fSide
886 fLowerEndcap->SetNeighbours( fSide180, fSide
887
888 }
889
890 //============================================
891 //* GetCubicVolume ---------------------------
892
893 G4double G4VTwistedFaceted::GetCubicVolume()
894 {
895 if(fCubicVolume == 0.)
896 {
897 fCubicVolume = ((fDx1 + fDx2 + fDx3 + fDx4
898 (fDx4 + fDx3 - fDx2 - fDx1
899 }
900 return fCubicVolume;
901 }
902
903 //============================================
904 //* GetLateralFaceArea -----------------------
905
906 G4double
907 G4VTwistedFaceted::GetLateralFaceArea(const G4
908 const G4
909 const G4
910 const G4
911 {
912 constexpr G4int NSTEP = 100;
913 constexpr G4double dt = 1./NSTEP;
914
915 G4double h = 2.*fDz;
916 G4double hh = h*h;
917 G4double hTanTheta = h*std::tan(fTheta);
918 G4double x1 = p1.x();
919 G4double y1 = p1.y();
920 G4double x21 = p2.x() - p1.x();
921 G4double y21 = p2.y() - p1.y();
922 G4double x31 = p3.x() - p1.x();
923 G4double y31 = p3.y() - p1.y();
924 G4double x42 = p4.x() - p2.x();
925 G4double y42 = p4.y() - p2.y();
926 G4double x43 = p4.x() - p3.x();
927 G4double y43 = p4.y() - p3.y();
928
929 // check if face is plane (just in case)
930 G4double lmax = std::max(std::max(std::abs(x
931 std::max(std::abs(x
932 G4double eps = lmax*kCarTolerance;
933 if (std::abs(fPhiTwist) < kCarTolerance &&
934 std::abs(x21*y43 - y21*x43) < eps)
935 {
936 G4double x0 = hTanTheta*std::cos(fPhi);
937 G4double y0 = hTanTheta*std::sin(fPhi);
938 G4ThreeVector A(p4.x() - p1.x() + x0, p4.y
939 G4ThreeVector B(p3.x() - p2.x() + x0, p3.y
940 return (A.cross(B)).mag()*0.5;
941 }
942
943 // twisted face
944 G4double area = 0.;
945 for (G4int i = 0; i < NSTEP; ++i)
946 {
947 G4double t = (i + 0.5)*dt;
948 G4double I = x21 + (x42 - x31)*t;
949 G4double J = y21 + (y42 - y31)*t;
950 G4double II = I*I;
951 G4double JJ = J*J;
952 G4double IIJJ = hh*(I*I + J*J);
953
954 G4double ang = fPhi + fPhiTwist*(0.5 - t);
955 G4double A = fPhiTwist*(II + JJ) + x21*y43
956 G4double B = fPhiTwist*(I*(x1 + x31*t) + J
957 hTanTheta*(I*std::sin(ang) - J*std::cos(
958 (I*y31 - J*x31);
959
960 G4double invAA = 1./(A*A);
961 G4double sqrtAA = 2.*std::abs(A);
962 G4double invSqrtAA = 1./sqrtAA;
963
964 G4double aa = A*A;
965 G4double bb = 2.*A*B;
966 G4double cc = IIJJ + B*B;
967
968 G4double R1 = std::sqrt(aa + bb + cc);
969 G4double R0 = std::sqrt(cc);
970 G4double log1 = std::log(std::abs(sqrtAA*R
971 G4double log0 = std::log(std::abs(sqrtAA*R
972
973 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) +
974 }
975 return area*dt;
976 }
977
978 //============================================
979 //* GetSurfaceArea ---------------------------
980
981 G4double G4VTwistedFaceted::GetSurfaceArea()
982 {
983 if (fSurfaceArea == 0.)
984 {
985 G4TwoVector vv[8];
986 vv[0] = G4TwoVector(-fDx1 - fDy1*fTAlph,-f
987 vv[1] = G4TwoVector( fDx1 - fDy1*fTAlph,-f
988 vv[2] = G4TwoVector(-fDx2 + fDy1*fTAlph, f
989 vv[3] = G4TwoVector( fDx2 + fDy1*fTAlph, f
990 vv[4] = G4TwoVector(-fDx3 - fDy2*fTAlph,-f
991 vv[5] = G4TwoVector( fDx3 - fDy2*fTAlph,-f
992 vv[6] = G4TwoVector(-fDx4 + fDy2*fTAlph, f
993 vv[7] = G4TwoVector( fDx4 + fDy2*fTAlph, f
994 fSurfaceArea = 2.*(fDy1*(fDx1 + fDx2) + fD
995 GetLateralFaceArea(vv[0], vv[1], vv[4],
996 GetLateralFaceArea(vv[1], vv[3], vv[5],
997 GetLateralFaceArea(vv[3], vv[2], vv[7],
998 GetLateralFaceArea(vv[2], vv[0], vv[6],
999 }
1000 return fSurfaceArea;
1001 }
1002
1003 //===========================================
1004 //* GetEntityType ---------------------------
1005
1006 G4GeometryType G4VTwistedFaceted::GetEntityTy
1007 {
1008 return {"G4VTwistedFaceted"};
1009 }
1010
1011
1012 //===========================================
1013 //* GetPolyhedron ---------------------------
1014
1015 G4Polyhedron* G4VTwistedFaceted::GetPolyhedro
1016 {
1017 if (fpPolyhedron == nullptr ||
1018 fRebuildPolyhedron ||
1019 fpPolyhedron->GetNumberOfRotationStepsA
1020 fpPolyhedron->GetNumberOfRotationSteps(
1021 {
1022 G4AutoLock l(&polyhedronMutex);
1023 delete fpPolyhedron;
1024 fpPolyhedron = CreatePolyhedron();
1025 fRebuildPolyhedron = false;
1026 l.unlock();
1027 }
1028
1029 return fpPolyhedron;
1030 }
1031
1032
1033 //===========================================
1034 //* GetPointInSolid -------------------------
1035
1036 G4ThreeVector G4VTwistedFaceted::GetPointInSo
1037 {
1038
1039
1040 // this routine is only used for a test
1041 // can be deleted ...
1042
1043 if ( z == fDz ) z -= 0.1*fDz ;
1044 if ( z == -fDz ) z += 0.1*fDz ;
1045
1046 G4double phi = z/(2*fDz)*fPhiTwist ;
1047
1048 return { fdeltaX * phi/fPhiTwist, fdeltaY *
1049 }
1050
1051
1052 //===========================================
1053 //* GetPointOnSurface -----------------------
1054
1055 G4ThreeVector G4VTwistedFaceted::GetPointOnSu
1056 {
1057
1058 G4double phi = G4RandFlat::shoot(-fPhiTwist
1059 G4double u , umin, umax ; // variable for
1060 G4double y ; // variable for
1061
1062 // Compute the areas. Attention: Only corre
1063 // where the twisting is done along the z-a
1064 // the computed surface area is more diffic
1065 // does not affect the tracking through the
1066
1067 G4double a1 = fSide0->GetSurfaceArea();
1068 G4double a2 = fSide90->GetSurfaceArea();
1069 G4double a3 = fSide180->GetSurfaceArea() ;
1070 G4double a4 = fSide270->GetSurfaceArea() ;
1071 G4double a5 = fLowerEndcap->GetSurfaceArea(
1072 G4double a6 = fUpperEndcap->GetSurfaceArea(
1073
1074 #ifdef G4TWISTDEBUG
1075 G4cout << "Surface 0 deg = " << a1 << G4e
1076 G4cout << "Surface 90 deg = " << a2 << G4e
1077 G4cout << "Surface 180 deg = " << a3 << G4e
1078 G4cout << "Surface 270 deg = " << a4 << G4e
1079 G4cout << "Surface Lower = " << a5 << G4e
1080 G4cout << "Surface Upper = " << a6 << G4e
1081 #endif
1082
1083 G4double chose = G4RandFlat::shoot(0.,a1 +
1084
1085 if(chose < a1)
1086 {
1087 umin = fSide0->GetBoundaryMin(phi) ;
1088 umax = fSide0->GetBoundaryMax(phi) ;
1089 u = G4RandFlat::shoot(umin,umax) ;
1090
1091 return fSide0->SurfacePoint(phi, u, true
1092 }
1093
1094 else if( (chose >= a1) && (chose < a1 + a2
1095 {
1096 umin = fSide90->GetBoundaryMin(phi) ;
1097 umax = fSide90->GetBoundaryMax(phi) ;
1098
1099 u = G4RandFlat::shoot(umin,umax) ;
1100
1101 return fSide90->SurfacePoint(phi, u, true
1102 }
1103 else if( (chose >= a1 + a2 ) && (chose < a1
1104 {
1105 umin = fSide180->GetBoundaryMin(phi) ;
1106 umax = fSide180->GetBoundaryMax(phi) ;
1107 u = G4RandFlat::shoot(umin,umax) ;
1108
1109 return fSide180->SurfacePoint(phi, u, tru
1110 }
1111 else if( (chose >= a1 + a2 + a3 ) && (chos
1112 {
1113 umin = fSide270->GetBoundaryMin(phi) ;
1114 umax = fSide270->GetBoundaryMax(phi) ;
1115 u = G4RandFlat::shoot(umin,umax) ;
1116 return fSide270->SurfacePoint(phi, u, tru
1117 }
1118 else if( (chose >= a1 + a2 + a3 + a4 ) &&
1119 {
1120 y = G4RandFlat::shoot(-fDy1,fDy1) ;
1121 umin = fLowerEndcap->GetBoundaryMin(y) ;
1122 umax = fLowerEndcap->GetBoundaryMax(y) ;
1123 u = G4RandFlat::shoot(umin,umax) ;
1124
1125 return fLowerEndcap->SurfacePoint(u,y,tru
1126 }
1127 else
1128 {
1129 y = G4RandFlat::shoot(-fDy2,fDy2) ;
1130 umin = fUpperEndcap->GetBoundaryMin(y) ;
1131 umax = fUpperEndcap->GetBoundaryMax(y) ;
1132 u = G4RandFlat::shoot(umin,umax) ;
1133
1134 return fUpperEndcap->SurfacePoint(u,y,tru
1135
1136 }
1137 }
1138
1139
1140 //===========================================
1141 //* CreatePolyhedron ------------------------
1142
1143 G4Polyhedron* G4VTwistedFaceted::CreatePolyhe
1144 {
1145 // number of meshes
1146 const G4int k =
1147 G4int(G4Polyhedron::GetNumberOfRotationStep
1148 std::abs(fPhiTwist) / twopi) + 2;
1149 const G4int n = k;
1150
1151 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k
1152 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1
1153
1154 auto ph = new G4Polyhedron;
1155 typedef G4double G4double3[3];
1156 typedef G4int G4int4[4];
1157 auto xyz = new G4double3[nnodes]; // numbe
1158 auto faces = new G4int4[nfaces] ; // numbe
1159
1160 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1161 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1162 fSide270->GetFacets(k,n,xyz,faces,2) ;
1163 fSide0->GetFacets(k,n,xyz,faces,3) ;
1164 fSide90->GetFacets(k,n,xyz,faces,4) ;
1165 fSide180->GetFacets(k,n,xyz,faces,5) ;
1166
1167 ph->createPolyhedron(nnodes,nfaces,xyz,face
1168
1169 return ph;
1170 }
1171