| 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 // G4TwistTubsHypeSide implementation
27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch),
30 // from original version in Jupi
31 // -------------------------------------------
32
33 #include "G4TwistTubsHypeSide.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4GeometryTolerance.hh"
36
37 //============================================
38 //* constructors -----------------------------
39
40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const
41 const
42 const
43 const
44 const
45 const
46 const
47 const
48 const
49
50
51
52
53 : G4VTwistSurface(name, rot, tlate, handedne
54 axis0min, axis1min, axis0m
55 fKappa(kappa), fTanStereo(tanstereo),
56 fTan2Stereo(tanstereo*tanstereo), fR0(r0),
57 {
58 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59 {
60 G4Exception("G4TwistTubsHypeSide::G4Twis
61 "GeomSolids0002", FatalError
62 "Should swap axis0 and axis1
63 }
64 fInside.gp.set(kInfinity, kInfinity, kInfin
65 fInside.inside = kOutside;
66 fIsValidNorm = false;
67 SetCorners();
68 SetBoundaries();
69 }
70
71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const
72
73
74
75
76
77
78
79
80
81
82
83 : G4VTwistSurface(name)
84 {
85
86 fHandedness = handedness; // +z = +ve, -z
87 fAxis[0] = kPhi;
88 fAxis[1] = kZAxis;
89 fAxisMin[0] = kInfinity; // we cann
90 fAxisMax[0] = kInfinity; // because
91 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1];
93 fKappa = Kappa;
94 fDPhi = DPhi ;
95
96 if (handedness < 0) // inner hyperbolic su
97 {
98 fTanStereo = TanInnerStereo;
99 fR0 = InnerRadius;
100 }
101 else // outer hyperbolic su
102 {
103 fTanStereo = TanOuterStereo;
104 fR0 = OuterRadius;
105 }
106 fTan2Stereo = fTanStereo * fTanStereo;
107 fR02 = fR0 * fR0;
108
109 fTrans.set(0, 0, 0);
110 fIsValidNorm = false;
111
112 fInside.gp.set(kInfinity, kInfinity, kInfin
113 fInside.inside = kOutside;
114
115 SetCorners(EndInnerRadius, EndOuterRadius,
116
117 SetBoundaries();
118 }
119
120 //============================================
121 //* Fake default constructor -----------------
122
123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __vo
124 : G4VTwistSurface(a)
125 {
126 }
127
128 //============================================
129 //* destructor -------------------------------
130
131 G4TwistTubsHypeSide::~G4TwistTubsHypeSide() =
132
133 //============================================
134 //* GetNormal --------------------------------
135
136 G4ThreeVector G4TwistTubsHypeSide::GetNormal(c
137
138 {
139 // GetNormal returns a normal vector at a s
140 // to surface) point at tmpxx.
141 // If isGlobal=true, it returns the normal
142
143 G4ThreeVector xx;
144 if (isGlobal)
145 {
146 xx = ComputeLocalPoint(tmpxx);
147 if ((xx - fCurrentNormal.p).mag() < 0.5
148 {
149 return ComputeGlobalDirection(fCurren
150 }
151 }
152 else
153 {
154 xx = tmpxx;
155 if (xx == fCurrentNormal.p)
156 {
157 return fCurrentNormal.normal;
158 }
159 }
160
161 fCurrentNormal.p = xx;
162
163 G4ThreeVector normal( xx.x(), xx.y(), -xx.z
164 normal *= fHandedness;
165 normal = normal.unit();
166
167 if (isGlobal)
168 {
169 fCurrentNormal.normal = ComputeLocalDire
170 }
171 else
172 {
173 fCurrentNormal.normal = normal;
174 }
175 return fCurrentNormal.normal;
176 }
177
178 //============================================
179 //* Inside() ---------------------------------
180
181 EInside G4TwistTubsHypeSide::Inside(const G4Th
182 {
183 // Inside returns
184 const G4double halftol
185 = 0.5 * G4GeometryTolerance::GetInstance(
186
187 if (fInside.gp == gp)
188 {
189 return fInside.inside;
190 }
191 fInside.gp = gp;
192
193 G4ThreeVector p = ComputeLocalPoint(gp);
194
195
196 if (p.mag() < DBL_MIN)
197 {
198 fInside.inside = kOutside;
199 return fInside.inside;
200 }
201
202 G4double rhohype = GetRhoAtPZ(p);
203 G4double distanceToOut = fHandedness * (rho
204 // +ve : inside
205
206 if (distanceToOut < -halftol)
207 {
208 fInside.inside = kOutside;
209 }
210 else
211 {
212 G4int areacode = GetAreaCode(p);
213 if (IsOutside(areacode))
214 {
215 fInside.inside = kOutside;
216 }
217 else if (IsBoundary(areacode))
218 {
219 fInside.inside = kSurface;
220 }
221 else if (IsInside(areacode))
222 {
223 if (distanceToOut <= halftol)
224 {
225 fInside.inside = kSurface;
226 }
227 else
228 {
229 fInside.inside = kInside;
230 }
231 }
232 else
233 {
234 G4cout << "WARNING - G4TwistTubsHypeS
235 << " Invalid option !
236 << " name, areacode,
237 << GetName() << ", " << std::h
238 << std::dec << ", " << distanc
239 }
240 }
241
242 return fInside.inside;
243 }
244
245 //============================================
246 //* DistanceToSurface ------------------------
247
248 G4int G4TwistTubsHypeSide::DistanceToSurface(c
249 c
250
251
252
253
254
255 {
256 // Decide if and where a line intersects wi
257 // surface (of infinite extent)
258 //
259 // Arguments:
260 // p - (in) Point on trajectory
261 // v - (in) Vector along trajecto
262 // r2 - (in) Square of radius at z
263 // tan2phi - (in) std::tan(stereo)**2
264 // s - (out) Up to two points of
265 // intersection point i
266 // two intersections, s
267 // Returns:
268 // The number of intersections. If 0, t
269 //
270 //
271 // Equation of a line:
272 //
273 // x = x0 + s*tx y = y0 + s*ty
274 //
275 // Equation of a hyperbolic surface:
276 //
277 // x**2 + y**2 = r**2 + (z*tanPhi)**2
278 //
279 // Solution is quadratic:
280 //
281 // a*s**2 + b*s + c = 0
282 //
283 // where:
284 //
285 // a = tx**2 + ty**2 - (tz*tanPhi)**2
286 //
287 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2
288 //
289 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)*
290 //
291
292 fCurStatWithV.ResetfDone(validate, &gp, &gv
293
294 if (fCurStatWithV.IsDone())
295 {
296 for (G4int i=0; i<fCurStatWithV.GetNXX()
297 {
298 gxx[i] = fCurStatWithV.GetXX(i);
299 distance[i] = fCurStatWithV.GetDistan
300 areacode[i] = fCurStatWithV.GetAreaco
301 isvalid[i] = fCurStatWithV.IsValid(i
302 }
303 return fCurStatWithV.GetNXX();
304 }
305 else // initialize
306 {
307 for (auto i=0; i<2; ++i)
308 {
309 distance[i] = kInfinity;
310 areacode[i] = sOutside;
311 isvalid[i] = false;
312 gxx[i].set(kInfinity, kInfinity, kInf
313 }
314 }
315
316 G4ThreeVector p = ComputeLocalPoint(gp);
317 G4ThreeVector v = ComputeLocalDirection(gv)
318 G4ThreeVector xx[2];
319
320 //
321 // special case! p is on origin.
322 //
323
324 if (p.mag() == 0)
325 {
326 // p is origin.
327 // unique solution of 2-dimension questi
328 // Equations:
329 // r^2 = fR02 + z^2*fTan2Stere0
330 // r = beta*z
331 // where
332 // beta = vrho / vz
333 // Solution (z value of intersection poi
334 // xxz = +- std::sqrt (fR02 / (beta^2
335 //
336
337 G4double vz = v.z();
338 G4double absvz = std::fabs(vz);
339 G4double vrho = v.getRho();
340 G4double vslope = vrho/vz;
341 G4double vslope2 = vslope * vslope;
342 if (vrho == 0 || (vrho/absvz) <= (absvz*
343 {
344 // vz/vrho is bigger than slope of as
345 distance[0] = kInfinity;
346 fCurStatWithV.SetCurrentStatus(0, gxx
347 isvali
348 return 0;
349 }
350
351 if (vz != 0.0)
352 {
353 G4double xxz = std::sqrt(fR02 / (vsl
354 * (vz / std::fabs(vz))
355 G4double t = xxz / vz;
356 xx[0].set(t*v.x(), t*v.y(), xxz);
357 }
358 else
359 {
360 // p.z = 0 && v.z =0
361 xx[0].set(v.x()*fR0, v.y()*fR0, 0);
362 }
363 distance[0] = xx[0].mag();
364 gxx[0] = ComputeGlobalPoint(xx[0]);
365
366 if (validate == kValidateWithTol)
367 {
368 areacode[0] = GetAreaCode(xx[0]);
369 if (!IsOutside(areacode[0]))
370 {
371 if (distance[0] >= 0) isvalid[0] =
372 }
373 }
374 else if (validate == kValidateWithoutTol
375 {
376 areacode[0] = GetAreaCode(xx[0], fals
377 if (IsInside(areacode[0]))
378 {
379 if (distance[0] >= 0) isvalid[0] =
380 }
381 }
382 else // kDontValidate
383 {
384 areacode[0] = sInside;
385 if (distance[0] >= 0) isvalid[0] =
386 }
387
388 fCurStatWithV.SetCurrentStatus(0, gxx[0]
389 isvalid[0
390 return 1;
391 }
392
393 //
394 // special case end.
395 //
396
397 G4double a = v.x()*v.x() + v.y()*v.y() - v.
398 G4double b = 2.0
399 * ( p.x() * v.x() + p.y() * v.y(
400 G4double c = p.x()*p.x() + p.y()*p.y() - fR
401 G4double D = b*b - 4*a*c; //discri
402 G4int vout = 0;
403
404 if (std::fabs(a) < DBL_MIN)
405 {
406 if (std::fabs(b) > DBL_MIN) /
407 {
408 distance[0] = -c/b;
409 xx[0] = p + distance[0]*v;
410 gxx[0] = ComputeGlobalPoint(xx[0]);
411
412 if (validate == kValidateWithTol)
413 {
414 areacode[0] = GetAreaCode(xx[0]);
415 if (!IsOutside(areacode[0]))
416 {
417 if (distance[0] >= 0) isvalid[0
418 }
419 }
420 else if (validate == kValidateWithout
421 {
422 areacode[0] = GetAreaCode(xx[0], f
423 if (IsInside(areacode[0]))
424 {
425 if (distance[0] >= 0) isvalid[0
426 }
427 }
428 else // kDontValidate
429 {
430 areacode[0] = sInside;
431 if (distance[0] >= 0) isvalid[0
432 }
433
434 fCurStatWithV.SetCurrentStatus(0, gxx
435 isvali
436 vout = 1;
437 }
438 else
439 {
440 // if a=b=0 and c != 0, p is origin an
441 // if a=b=c=0, p is on surface and v i
442 // return distance = infinity
443
444 fCurStatWithV.SetCurrentStatus(0, gxx
445 isvali
446 vout = 0;
447 }
448 }
449 else if (D > DBL_MIN) // double so
450 {
451 D = std::sqrt(D);
452 G4double factor = 0.5/a;
453 G4double tmpdist[2] = {kInfinity, k
454 G4ThreeVector tmpxx[2] ;
455 G4int tmpareacode[2] = {sOutside
456 G4bool tmpisvalid[2] = {false, f
457
458 for (auto i=0; i<2; ++i)
459 {
460 tmpdist[i] = factor*(-b - D);
461 D = -D;
462 tmpxx[i] = p + tmpdist[i]*v;
463
464 if (validate == kValidateWithTol)
465 {
466 tmpareacode[i] = GetAreaCode(tmpxx
467 if (!IsOutside(tmpareacode[i]))
468 {
469 if (tmpdist[i] >= 0) tmpisvalid
470 continue;
471 }
472 }
473 else if (validate == kValidateWithout
474 {
475 tmpareacode[i] = GetAreaCode(tmpxx
476 if (IsInside(tmpareacode[i]))
477 {
478 if (tmpdist[i] >= 0) tmpisvalid
479 continue;
480 }
481 }
482 else // kDontValidate
483 {
484 tmpareacode[i] = sInside;
485 if (tmpdist[i] >= 0) tmpisvalid
486 continue;
487 }
488 }
489
490 if (tmpdist[0] <= tmpdist[1])
491 {
492 distance[0] = tmpdist[0];
493 distance[1] = tmpdist[1];
494 xx[0] = tmpxx[0];
495 xx[1] = tmpxx[1];
496 gxx[0] = ComputeGlobalPoint(tmp
497 gxx[1] = ComputeGlobalPoint(tmp
498 areacode[0] = tmpareacode[0];
499 areacode[1] = tmpareacode[1];
500 isvalid[0] = tmpisvalid[0];
501 isvalid[1] = tmpisvalid[1];
502 }
503 else
504 {
505 distance[0] = tmpdist[1];
506 distance[1] = tmpdist[0];
507 xx[0] = tmpxx[1];
508 xx[1] = tmpxx[0];
509 gxx[0] = ComputeGlobalPoint(tmp
510 gxx[1] = ComputeGlobalPoint(tmp
511 areacode[0] = tmpareacode[1];
512 areacode[1] = tmpareacode[0];
513 isvalid[0] = tmpisvalid[1];
514 isvalid[1] = tmpisvalid[0];
515 }
516
517 fCurStatWithV.SetCurrentStatus(0, gxx[0]
518 isvalid[0
519 fCurStatWithV.SetCurrentStatus(1, gxx[1]
520 isvalid[1
521 vout = 2;
522 }
523 else
524 {
525 // if D<0, no solution
526 // if D=0, just grazing the surfaces, re
527
528 fCurStatWithV.SetCurrentStatus(0, gxx[0]
529 isvalid[0
530 vout = 0;
531 }
532 return vout;
533 }
534
535 //============================================
536 //* DistanceToSurface ------------------------
537
538 G4int G4TwistTubsHypeSide::DistanceToSurface(c
539
540
541
542 {
543 // Find the approximate distance of a poin
544 // The distance must be an underestimate.
545 // It will also be nice (although not nece
546 // always finite no matter how close the p
547 //
548 // We arranged G4Hype::ApproxDistOutside a
549 // for this function. See these discriptio
550
551 const G4double halftol
552 = 0.5 * G4GeometryTolerance::GetInstance(
553
554 fCurStat.ResetfDone(kDontValidate, &gp);
555
556 if (fCurStat.IsDone())
557 {
558 for (G4int i=0; i<fCurStat.GetNXX(); ++i
559 {
560 gxx[i] = fCurStat.GetXX(i);
561 distance[i] = fCurStat.GetDistance(i)
562 areacode[i] = fCurStat.GetAreacode(i)
563 }
564 return fCurStat.GetNXX();
565 }
566 else // initialize
567 {
568 for (auto i=0; i<2; ++i)
569 {
570 distance[i] = kInfinity;
571 areacode[i] = sOutside;
572 gxx[i].set(kInfinity, kInfinity, kInf
573 }
574 }
575
576
577 G4ThreeVector p = ComputeLocalPoint(gp);
578 G4ThreeVector xx;
579
580 //
581 // special case!
582 // If p is on surface, return distance = 0
583 //
584 G4ThreeVector lastgxx[2];
585 for (auto i=0; i<2; ++i)
586 {
587 lastgxx[i] = fCurStatWithV.GetXX(i);
588 }
589
590 if ((gp - lastgxx[0]).mag() < halftol || (g
591 {
592 // last winner, or last poststep point i
593 xx = p;
594 gxx[0] = gp;
595 distance[0] = 0;
596
597 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], dis
599 isvalid, 1, kD
600
601 return 1;
602 }
603 //
604 // special case end
605 //
606
607 G4double prho = p.getRho();
608 G4double pz = std::fabs(p.z());
609 G4double r1 = std::sqrt(fR02 + pz *
610
611 G4ThreeVector pabsz(p.x(), p.y(), pz);
612
613 if (prho > r1 + halftol) // p is outside
614 {
615 // First point xx1
616 G4double t = r1 / prho;
617 G4ThreeVector xx1(t * pabsz.x(), t * pab
618
619 // Second point xx2
620 G4double z2 = (prho * fTanStereo + pz) /
621 G4double r2 = std::sqrt(fR02 + z2 * z2 *
622 t = r2 / prho;
623 G4ThreeVector xx2(t * pabsz.x(), t * pab
624
625 G4double len = (xx2 - xx1).mag();
626 if (len < DBL_MIN)
627 {
628 // xx2 = xx1?? I guess we
629 // must have really bracketed the nor
630 distance[0] = (pabsz - xx1).mag();
631 xx = xx1;
632 }
633 else
634 {
635 distance[0] = DistanceToLine(pabsz, x
636 }
637
638 }
639 else if (prho < r1 - halftol) // p is insi
640 {
641 // First point xx1
642 G4double t;
643 G4ThreeVector xx1;
644 if (prho < DBL_MIN)
645 {
646 xx1.set(r1, 0. , pz);
647 }
648 else
649 {
650 t = r1 / prho;
651 xx1.set(t * pabsz.x(), t * pabsz.y()
652 }
653
654 // dr, dz is tangential vector of Hyparb
655 // dr = r, dz = z*tan2stereo
656 G4double dr = pz * fTan2Stereo;
657 G4double dz = r1;
658 G4double tanbeta = dr / dz;
659 G4double pztanbeta = pz * tanbeta;
660
661 // Second point xx2
662 // xx2 is intersection between x-axis an
663 G4double r2 = r1 - pztanbeta;
664 G4ThreeVector xx2;
665 if (prho < DBL_MIN)
666 {
667 xx2.set(r2, 0. , 0.);
668 }
669 else
670 {
671 t = r2 / prho;
672 xx2.set(t * pabsz.x(), t * pabsz.y()
673 }
674
675 G4ThreeVector d = xx2 - xx1;
676 distance[0] = DistanceToLine(pabsz, xx1,
677
678 }
679 else // p is on Hyperbolic surface.
680 {
681 distance[0] = 0;
682 xx.set(p.x(), p.y(), pz);
683 }
684
685 if (p.z() < 0)
686 {
687 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.
688 xx = tmpxx;
689 }
690
691 gxx[0] = ComputeGlobalPoint(xx);
692 areacode[0] = sInside;
693 G4bool isvalid = true;
694 fCurStat.SetCurrentStatus(0, gxx[0], distan
695 isvalid, 1, kDont
696 return 1;
697 }
698
699 //============================================
700 //* GetAreaCode ------------------------------
701
702 G4int G4TwistTubsHypeSide::GetAreaCode(const G
703 G
704 {
705 const G4double ctol = 0.5 * kCarTolerance;
706 G4int areacode = sInside;
707
708 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis
709 {
710 G4int zaxis = 1;
711
712 if (withTol)
713 {
714 G4bool isoutside = false;
715 G4int phiareacode = GetAreaCodeIn
716 G4bool isoutsideinphi = IsOutside(phi
717
718 // test boundary of phiaxis
719
720 if ((phiareacode & sAxisMin) == sAxis
721 {
722 areacode |= (sAxis0 & (sAxisPhi |
723 if (isoutsideinphi) isoutside = tr
724 }
725 else if ((phiareacode & sAxisMax) ==
726 {
727 areacode |= (sAxis0 & (sAxisPhi |
728 if (isoutsideinphi) isoutside = tr
729 }
730
731 // test boundary of zaxis
732
733 if (xx.z() < fAxisMin[zaxis] + ctol)
734 {
735 areacode |= (sAxis1 & (sAxisZ | sA
736 if ((areacode & sBoundary) != 0)
737 else areaco
738
739 if (xx.z() <= fAxisMin[zaxis] - ct
740
741 }
742 else if (xx.z() > fAxisMax[zaxis] - c
743 {
744 areacode |= (sAxis1 & (sAxisZ | sA
745 if ((areacode & sBoundary) != 0)
746 else areaco
747
748 if (xx.z() >= fAxisMax[zaxis] + ct
749 }
750
751 // if isoutside = true, clear sInside
752 // if not on boundary, add boundary i
753
754 if (isoutside)
755 {
756 G4int tmpareacode = areacode & (~s
757 areacode = tmpareacode;
758 }
759 else if ((areacode & sBoundary) != sB
760 {
761 areacode |= (sAxis0 & sAxisPhi) |
762 }
763
764 return areacode;
765 }
766 else
767 {
768 G4int phiareacode = GetAreaCodeInPhi(
769
770 // test boundary of z-axis
771
772 if (xx.z() < fAxisMin[zaxis])
773 {
774 areacode |= (sAxis1 & (sAxisZ | sA
775
776 }
777 else if (xx.z() > fAxisMax[zaxis])
778 {
779 areacode |= (sAxis1 & (sAxisZ | sA
780 }
781
782 // boundary of phi-axis
783
784 if (phiareacode == sAxisMin)
785 {
786 areacode |= (sAxis0 & (sAxisPhi |
787 if ((areacode & sBoundary) != 0)
788 else areaco
789
790 }
791 else if (phiareacode == sAxisMax)
792 {
793 areacode |= (sAxis0 & (sAxisPhi |
794 if ((areacode & sBoundary) != 0)
795 else areaco
796 }
797
798 // if not on boundary, add boundary i
799
800 if ((areacode & sBoundary) != sBounda
801 {
802 areacode |= (sAxis0 & sAxisPhi) |
803 }
804 return areacode;
805 }
806 }
807 else
808 {
809 std::ostringstream message;
810 message << "Feature NOT implemented !" <
811 << " fAxis[0] = " << fAxi
812 << " fAxis[1] = " << fAxi
813 G4Exception("G4TwistTubsHypeSide::GetAre
814 "GeomSolids0001", FatalExcep
815 }
816 return areacode;
817 }
818
819 //============================================
820 //* GetAreaCodeInPhi -------------------------
821
822 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(co
823
824 {
825
826 G4ThreeVector lowerlimit; // lower phi-boun
827 G4ThreeVector upperlimit; // upper phi-boun
828 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxis
829 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxis
830
831 G4int areacode = sInside;
832 G4bool isoutside = false;
833
834 if (withTol)
835 {
836 if (AmIOnLeftSide(xx, lowerlimit) >= 0)
837 {
838 areacode |= (sAxisMin | sBoundary);
839 if (AmIOnLeftSide(xx, lowerlimit) > 0
840
841 }
842 else if (AmIOnLeftSide(xx, upperlimit) <
843 {
844 areacode |= (sAxisMax | sBoundary);
845 if (AmIOnLeftSide(xx, upperlimit) < 0
846 }
847
848 // if isoutside = true, clear inside bit
849
850 if (isoutside)
851 {
852 G4int tmpareacode = areacode & (~sIns
853 areacode = tmpareacode;
854 }
855 }
856 else
857 {
858 if (AmIOnLeftSide(xx, lowerlimit, false)
859 {
860 areacode |= (sAxisMin | sBoundary);
861 }
862 else if (AmIOnLeftSide(xx, upperlimit, f
863 {
864 areacode |= (sAxisMax | sBoundary);
865 }
866 }
867
868 return areacode;
869 }
870
871 //============================================
872 //* SetCorners(EndInnerRadius, EndOuterRadius,
873
874 void G4TwistTubsHypeSide::SetCorners( G4double
875 G4double
876 G4double
877 G4double
878 G4double
879 {
880 // Set Corner points in local coodinate.
881
882 if (fAxis[0] == kPhi && fAxis[1] == kZAxis)
883
884 G4double endRad[2];
885 G4double halfdphi = 0.5*DPhi;
886
887 for (auto i=0; i<2; ++i) // i=0,1 : -ve
888 {
889 endRad[i] = (fHandedness == 1 ? EndOut
890 }
891
892 G4int zmin = 0 ; // at -ve z
893 G4int zmax = 1 ; // at +ve z
894
895 G4double x, y, z;
896
897 // corner of Axis0min and Axis1min
898 x = endRad[zmin]*std::cos(endPhi[zmin] -
899 y = endRad[zmin]*std::sin(endPhi[zmin] -
900 z = endZ[zmin];
901 SetCorner(sC0Min1Min, x, y, z);
902
903 // corner of Axis0max and Axis1min
904 x = endRad[zmin]*std::cos(endPhi[zmin] +
905 y = endRad[zmin]*std::sin(endPhi[zmin] +
906 z = endZ[zmin];
907 SetCorner(sC0Max1Min, x, y, z);
908
909 // corner of Axis0max and Axis1max
910 x = endRad[zmax]*std::cos(endPhi[zmax] +
911 y = endRad[zmax]*std::sin(endPhi[zmax] +
912 z = endZ[zmax];
913 SetCorner(sC0Max1Max, x, y, z);
914
915 // corner of Axis0min and Axis1max
916 x = endRad[zmax]*std::cos(endPhi[zmax] -
917 y = endRad[zmax]*std::sin(endPhi[zmax] -
918 z = endZ[zmax];
919 SetCorner(sC0Min1Max, x, y, z);
920
921 }
922 else
923 {
924 std::ostringstream message;
925 message << "Feature NOT implemented !" <
926 << " fAxis[0] = " << fAxi
927 << " fAxis[1] = " << fAxi
928 G4Exception("G4TwistTubsHypeSide::SetCor
929 "GeomSolids0001", FatalExcep
930 }
931 }
932
933 //============================================
934 //* SetCorners() -----------------------------
935
936 void G4TwistTubsHypeSide::SetCorners()
937 {
938 G4Exception("G4TwistTubsHypeSide::SetCorner
939 "GeomSolids0001", FatalExceptio
940 "Method NOT implemented !");
941 }
942
943 //============================================
944 //* SetBoundaries() --------------------------
945
946 void G4TwistTubsHypeSide::SetBoundaries()
947 {
948 // Set direction-unit vector of phi-boundar
949 // sAxis0 must be kPhi.
950 // This fanction set lower phi-boundary and
951
952 if (fAxis[0] == kPhi && fAxis[1] == kZAxis)
953 {
954 G4ThreeVector direction;
955 // sAxis0 & sAxisMin
956 direction = GetCorner(sC0Min1Max) - GetC
957 direction = direction.unit();
958 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi
959 GetCorner(sC0Min1Min), sAxi
960
961 // sAxis0 & sAxisMax
962 direction = GetCorner(sC0Max1Max) - GetC
963 direction = direction.unit();
964 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa
965 GetCorner(sC0Max1Min), sAxis
966
967 // sAxis1 & sAxisMin
968 direction = GetCorner(sC0Max1Min) - GetC
969 direction = direction.unit();
970 SetBoundary(sAxis1 & (sAxisZ | sAxisMin)
971 GetCorner(sC0Min1Min), sAxi
972
973 // sAxis1 & sAxisMax
974 direction = GetCorner(sC0Max1Max) - GetC
975 direction = direction.unit();
976 SetBoundary(sAxis1 & (sAxisZ | sAxisMax)
977 GetCorner(sC0Min1Max), sAxis
978 }
979 else
980 {
981 std::ostringstream message;
982 message << "Feature NOT implemented !" <
983 << " fAxis[0] = " << fAxi
984 << " fAxis[1] = " << fAxi
985 G4Exception("G4TwistTubsHypeSide::SetBou
986 "GeomSolids0001", FatalExcep
987 }
988 }
989
990 //============================================
991 //* GetFacets() ------------------------------
992
993 void G4TwistTubsHypeSide::GetFacets( G4int k,
994 G4int fac
995 {
996 G4double z ; // the two parameters for t
997 G4double x,xmin,xmax ;
998
999 G4ThreeVector p ; // a point on the surface
1000
1001 G4int nnode ;
1002 G4int nface ;
1003
1004 // calculate the (n-1)*(k-1) vertices
1005
1006 for ( G4int i = 0 ; i<n ; ++i )
1007 {
1008 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin
1009
1010 for ( G4int j = 0 ; j<k ; ++j )
1011 {
1012 nnode = GetNode(i,j,k,n,iside) ;
1013
1014 xmin = GetBoundaryMin(z) ;
1015 xmax = GetBoundaryMax(z) ;
1016
1017 if (fHandedness < 0) // inner hyperbol
1018 {
1019 x = xmin + j*(xmax-xmin)/(k-1) ;
1020 }
1021 else // outer hyperbol
1022 {
1023 x = xmax - j*(xmax-xmin)/(k-1) ;
1024 }
1025
1026 p = SurfacePoint(x,z,true) ; // surfac
1027
1028 xyz[nnode][0] = p.x() ;
1029 xyz[nnode][1] = p.y() ;
1030 xyz[nnode][2] = p.z() ;
1031
1032 if ( i<n-1 && j<k-1 ) // clock wise
1033 {
1034 nface = GetFace(i,j,k,n,iside) ;
1035 faces[nface][0] = GetEdgeVisibility(i,j,k,n
1036 * ( GetNode(i ,j ,k
1037 faces[nface][1] = GetEdgeVisibility(i,j,k,n
1038 * ( GetNode(i+1,j ,k
1039 faces[nface][2] = GetEdgeVisibility(i,j,k,n
1040 * ( GetNode(i+1,j+1,k
1041 faces[nface][3] = GetEdgeVisibility(i,j,k,n
1042 * ( GetNode(i ,j+1,k
1043 }
1044 }
1045 }
1046 }
1047